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ABSTRACT 

We studied roles of a turbulent resistivity in the core-collapse of a strongly magnetized massive star, 
carrying out 2D-resistive-MHD simulations. The three cases with different initial strengths of magnetic 
field and rotation are investigated; 1. strongly magnetized rotating core; 2. moderately magnetized 
rotating core; 3. very strongly magnetized non-rotating core. In each case, both an ideal-MHD model 
and resistive-MHD models are computed. As a result of computations, each model shows a matter 
eruption helped by a magnetic acceleration (and also by a centrifugal acceleration in the rotating 
cases). We found that a resistivity attenuates the explosion in case 1 and 2, while it enhances the 
explosion in case 3. We also found that in the rotating cases, main mechanisms for the amplification 
of a magnetic field in the post-bounce phase are an outward advection of magnetic field and a winding 
of poloidal magnetic field-lines by differential rotation, which are somewhat dampened down with the 
presence of a resistivity. Although the magnetorotational instability seems to occur in the rotating 
models, it will play only a minor role in a magnetic field amplification. Another impact of resistivity 
is that on the aspect ratio. In the rotating cases, a large aspect ratio of the ejected matters, > 2.5, 
attained in a ideal-MHD model is reduced to some extent in a resistive model. These results indicate 
that a resistivity possibly plays an important role in the dynamics of strongly magnetized supernovae. 
Subject headings: magnetohydrodynamics (MHD) — methods: numerical — stars: magnetars — 
supernovae: general 



1. INTRODUCTION 

Studies on magnetized core-collapse supernovae (CC- 
SNe) has gathered stream for past several years. Numer- 
ical simulations so far have shown that the presence of 
strong magnetic field together with rapid rotation results 
in a vigorous matter eruption accompanied by bipolar 
jets with an explosion energy of ~ 10 51 erg (magnetoro- 
tational explosion: LeBlanc & Wilson 1970; Symbalisty 
1984; Yamada & Sawai 2004). One of the main driving 
forces is a toroidal magnetic pressure amplified by dif- 
ferential rotation that becomes intense after the collapse 
in the vicinity of proto-neutron star surface. This mech- 
anism requires a magnetic field strength of > 10 15 G 
after collapse, which is comparable to an inferred sur- 
face magnetic field of magnet ar candidates, soft-gamma 
repeaters (SGR) and anomalous X-ray pulsars (AXP). 
Magnetically-driven explosions may be related to such 
supernovae that produce magnetars. 

Most numerical simulations of magnetized core- 
collapse so far have been done in the regime of ideal 
magnetohydrodynamics (MHD) (e.g. Kotake et al. 2004; 
Sawai et al. 2005; Moiseenko et al. 2006; Burrows et 
al. 2007; Sawai et al. 2008; Takiwaki et al. 2009; Ober- 
gaulinger & Janka 2011; Endeve et al. 2012). The only 
exception is the numerical study by Guilet et al. (2011), 
in which a resistivity is introduced to 1D-MHD simu- 
lations investigating the dynamics of an Alfven surface 
in the context of core-collapse supernovae. Although a 
numerical computation with a finite differential scheme 
inevitably involves numerical diffusion, effects of electric 
resistivity on the dynamics have not been investigated 
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systematically. The reason why a resistivity has been ne- 
glected is that it would be infinitesimally small assuming 
that the source is the Coulomb scattering (Spitzer resis- 
tivity; Spitzer (1956)). The magnetic Reynolds num- 
ber, the ratio of the resistive timescale to dynamical 
timescale, estimated with a typical parameters of proto- 
neutron star surface is 
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where Z, T, L, and v are, respectively, an atomic charge 
number, temperature, scale length of magnetic field, and 
flow velocity. Since the magnetic Reynolds number is 
quite larger than unity, a resistivity apparently seems 
not important for the dynamics. 

However, it is uncertain whether the Coulomb scat- 
tering is the unique origin of resistivity in a supernova 
core, and there may exist other sources that give rise to a 
dynamically important value of resistivity. One of such 
candidates is a turbulence. In the collapsed core of a 
massive star, a convection, which occurs due to negative 
gradient of entropy and/or electron fraction, may play 
a key role to produce a turbulent state and a turbulent 
resistivity (and viscosity) along with that. Thompson et 
al. (2005) roughly estimated the amplitude of a turbulent 
viscosity arising from convective motions in a supernova 
core by the product of the correlation length and con- 
vective velocity, £ con ~ ^con/3. They found that £ con is 
around 10 13 — 10 14 cm 2 s -1 . Insofar as this level of esti- 
mation, the amplitude of a turbulent resistivity may be 
evaluated by the same formula and may be comparable 
to a turbulent viscosity, since magnetic field would have 
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a similar timescale and length-scale with those of veloc- 
ity in the present situation. Yoshizawa (1990) showed in 
the frame work of so-called two-scale direct-interaction 
approximation that the magnitude of a turbulent viscos- 
ity and turbulent resistivity are same order, albeit in the 
context of an incompressible MHD turbulence. With the 
above amplitude of the turbulent resistivity, the mag- 
netic Reynolds number becomes ~ 1 — 10 around the 
surface of a proto-neutron star, and then a resistivity is 
possibly important for the dynamics. 

In this paper, we investigate how a (turbulent) resis- 
tivity alters the dynamics of a magnetized core-collapse, 
paying particular attention to the explosion energy, the 
magnetic field amplification, and the aspect ratio of 
ejected matters. To this end we carried out axisymmet- 
ric 2D-resistive-MHD simulations of the core-collapse of 
a massive star, assuming a strong magnetic field and 
a large resistivity. A constant resistivity of 10 13 and 
10 14 cm 2 s _1 are taken according to the above discussion. 
Both rapidly rotating and non-rotating cores are studied. 
In the computations, we omitted any treatments of neu- 
trinos, and adopted a nuclear equation of state (EOS) 
produced by Shen et al. (1998a, b). 

Before proceeding to the next section, we go into a little 
more detail on a convection as a source of a turbulent 
resistivity, and also mention our position in choosing the 
initial strength of a magnetic field and a rotation. 

In a collapsed stellar core, there are mainly two con- 
vectively unstable regions (see e.g. Herant et al. (1994)). 
One is a region behind the shock surface, where a nega- 
tive entropy gradient is created as the shock surface prop- 
agates with its amplitude decreasing, and is maintained 
by a neutrino heating. The other is a region around the 
proto-neutron star surface, where a negative gradient of 
lepton fraction is created because a neutrino is easier to 
escape from the core for a lager radius. Since we do 
not deal with neutrinos, convections related to them are 
not captured in the computations. Only convection that 
seems to appear in our computation is one due to the 
shock propagation with a decreasing amplitude 4 . Never- 
theless, we assume all of the above convections as sources 
of a turbulent resistivity adopted in this study, since they 
will occur in the nature. What we have done here is to 
effectively introduce an impact of these convections to 
the simulations. In so doing, it does not seem quite im- 
portant whether they are properly captured in the simu- 
lations. Note that we do not consider that the standing 
accretion shock instability (SASI), which do not present 
in our computations, is one of the origins of a turbulent 
resistivity, because all of our models explode before this 
instability can develop (several 100 ms after bounce), ow- 
ing to strong magnetic fields initially assumed. 

Although a turbulent resistivity will appear only 
around the convectively unstable regions, in our compu- 
tations a constant resistivity is assumed everywhere ex- 
cept in the vicinity of the center (see § 3). Also, we should 
note that the estimation made by Thompson et al. (2005) 
is very uncertain. Moreover, strong magnetic fields ini- 
tially assumed may decrease the strength of a turbulent 

4 In each computation, we found that the square of Brunt- 
Baisala frequency is negative due to a negative entropy gradient 
in some locations behind the shock surface, and that a relatively 
large vorticity develops around there. 



resistivity. Therefore, the strengths of an adopted turbu- 
lent resistivity are perhaps too large to be realistic. How- 
ever, at present a probable value of a turbulent resistivity 
in a collapsed-stellar core is very unclear. Under such a 
circumstance, it is meaningful to parametrically study 
its effect with some possible values, and to grasp dynam- 
ical trends. We consider that the adopted strengths of a 
turbulent resistivity may be maximum possible values. 

A strongly magnetized core prior to collapse assumed 
in the present study is based on so-called "fossil-field hy- 
pothesis" which supposes that the progenitor of a mag- 
netar already has a magnetar-class magnetic flux dur- 
ing the main sequence stage. Assuming this hypothesis, 
Ferrario & Wickramasinghe (2006) have done popula- 
tion synthesis calculations from main sequence stars to 
neutron stars, to fit observational data of radio pulsars. 
Their calculation produces consistent number of magne- 
tars with those observed in the Galaxy, where both the 
age and rotational period of magnetars are taken into 
account. The fossil field hypothesis is also supported by 
observations. There are several O-type stars whose sur- 
face magnetic flux is inferred to be a magnetar-class; e.g. 
HD148937 (Wade et al. 2012) and HD19612 (Donati et 
al. 2006). Auriere et al. (2010) measured surface mag- 
netic fields of Betelgeuse, a red supergiant star, to be 
~ 1 G, which indicates a magnetar-class magnetic flux. 
Note, however, that the origin of a strong magnetic fields 
in magnetars is still controversial. Alternatively, they 
may be produced during the core-collapse by a dynamo 
mechanism (Thompson & Duncan 2003). 

Heger et al. (2005) found that the so-called Tayler- 
Spruit dynamo (Spruit 2002) drastically slows down the 
rotation of a star especially during an early phase of red 
super giant, where an angular momentum of the rapidly 
rotating helium core is transported into the slowly rotat- 
ing hydrogen envelope. According to their computation, 
an inferred rotational period of a pulser is ~ 10 ms for a 
15 Mq progenitor, in which the available rotational en- 
ergy is insufficient for the explosion. On the other hand, 
Woosley & Heger (2006) carried out stellar evolution 
computations of inherently rapid rotators, and showed 
that the rotation of a pre-supernova core could be fast. 
Due to the fast rotation, matters are almost completely 
mixed, and instead of forming a red supergiant it be- 
comes a compact helium core star, where a magnetic 
torque works less efficiently. One of their 16M mag- 
netic star models with the solar metalicity results in the 
expected pulsar rotation period of 2.3 ms, sufficient for a 
magnetorotational explosion. Note that the both works 
involve uncertainties about such as mass loss rate and 
multi-dimensional effects. At present, it is unclear either 
of slow or fast rotation is appropriate for the progenitor 
of a magnetar. Hence, in this study both rapidly rotating 
models and non-rotating models are investigated, where 
the latter, in effect, corresponds to a slow rotation case. 

The rest of this paper is organized as follows. We 
describe the governing equations and essentials of our 
resistive-MHD code in § 2, and computational setups in 
§ 3. The results are presented in § 4. Discussion and 
conclusion are given in § 5. 

2. GOVERNING EQUATIONS AND NUMERICAL SCHEMES 

In order to follow the dynamics of magnetized core- 
collapses with resistivity, the resistive-MHD equations 
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in which notations of the physical variables follow cus- 
tom. 

To solve above equations we have developed 2D- 
resistive-MHD code, " Yamazakura." This is a time ex- 
plicit, Eulerian code based on the high resolution central 
scheme formulated by Kurganov & Tadmor (2000) (Here 
after KT scheme). Below we briefly describe the features 
of Yamazakura. For sake of simplicity, we deal in the case 
where the equations are written in Cartesian coordinate 
with plane symmetry in z-direction. 

The KT scheme adopts a finite volume method to solve 
conservation equations. Although the induction equa- 
tions (4) apparently seem written in non-conservation 
forms, they are rewritten into conservation forms (Ziegler 
2004); 
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Then evolutionary Eqs. (1) — (3) and (7) are all written 
in conservation forms with source terms; 
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where a expression such as V • (f x , f y ) means df x /dx + 
df y /dy. The vectors in Eq. (8), each of which has eight 
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The KT scheme is written as follows: 
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where Axi = x i + 1/2 -x i _ 1/2 , Ayj = y j+ i/ 2 -yj-i/2, and 
an integer and half-integer subscript respectively means 
that a variable is evaluated at the numerical cell center 
and interface. The numerical fluxes F and G are for the 
non-resistive terms and resistive terms, respectively. The 
numerical fluxes of non-resistive terms are given by 
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The numerical fluxes of the resistive terms, G x and G y , 
are to appear later. 

In the original KT scheme, the interface values in 
Eq. (10), which have a superscript "+" or "— ", are eval- 
uated by a interpolation of conservative variables u: 
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Alternatively, we may be able to use interpolated 
interface values of the primitive variables q = 

(p, v x , v y , v z , e, B x , B y , B z ), 
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to evaluate the interface values of the conservative vari- 
ables in Eqs. (10), e.g., 

In Yamazakura, we adopt the latter prescription. For the 
calculation of the numerical derivatives, q and q , in 
Eqs. (12), we employ a minmod-like limiter suggested by 
Kurganov & Tadmor (2000). Note that the interpola- 
tions (12) are only used for the hydrodynamical quan- 
tities and the z-component of magnetic field, which are 
placed at the center of a numerical cell. Those for the x- 
component and ^/-component of magnetic field are given 
later. 

In Eqs. (10), a i + 1 / 2 j(t) and a^ J+1 / 2 (£) are the maxi- 
mum characteristic speeds, i.e. the sum of the fluid ve- 
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Figure 1. Positional relation of magnetic field vectors and nu- 
merical fluxes (or electric field vectors), E. 



locity and fast magnetosonic speed, which we evaluate 
using the interpolated primitive variables in Eqs. (12): 

a i+ i /2J (£)=max{a (g+ +1/2J (£)) ,a (V+i/2,j (*)) } , 

a*-i/2,iW=max{a (<?+i/2,j K>) > a (Ci/2jW) } ■ 

(13) 

In KT scheme, we only need to know these maximum 
characteristic speeds instead of carrying out a compli- 
cated characteristic decomposition for wave propaga- 
tions. 

The original KT scheme is formulated for uniform spa- 
tial cells. We have followed the procedure to deduce the 
KT scheme described in Kurganov & Tadmor (2000) with 
non- uniform cells, and found that the final semi-discrete 
form (9) — (11) is unchanged, except that the subscripts 
to Ax and Ay appear. 

In order to obtain the time evolution of conservative 
variables Uij(t), the semi-discrete equation (9) is time 
integrated utilizing a third order Runge-Kutta method 
according to Kurganov & Tadmor (2000). With this and 
the spatial interpolations in Eqs. (11), the original KT 
scheme is a third order in time and second order in space. 
In Appendix ("Linear Wave Propagation"), it is shown 
that Yamazakura performs, approximately, at least sec- 
ond order in time and second order in space, even though 
we adopt a non-uniform cell distribution and the spatial 
interpolations of the primitive variables (Eqs. (12)). 

In solving MHD equations, it is necessary to satisfy the 
divergence-free constraint of magnetic field. To accom- 
plish this, we apply a constraint transport (CT) method 
to KT scheme based on Ziegler (2004), extending it into 
resistive-MHD case. In a 3D-CT method, a magnetic 
field vector is placed at the center of a cubic cell inter- 
face while a numerical flux (or an electric field vector) 
is at a cell edge so that V • B does not evolve (Evans 
& Hawley 1988). The positional relation between them 
in 2D case are shown in Fig 1. Due to this placement, 
the semi-discrete equations and the numerical fluxes of 
the induction equations should be different from Eqs. (9) 
and (10). With numerical fluxes, E = F + G, the semi- 
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discrete form of the induction equations are written as 
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where F^ x,y ^ m ^ denotes the m-th component of the vec- 
tor F^ x,y \ The interpolations for the ^-component and 
^/-component of a magnetic field along the x-direction are 
given by 
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Those along the y-direction are given similar way to the 
above. Note again that interpolations of B z are same as 
Eq. (12), since in 2D it is defined at a cell center. 

In order to obtain the numerical fluxes of the resistive 
terms that appear in the energy equation (3) and the 
induction equations (7), an evaluation of current density 
is required, which we simply give by 
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By virtue of these definitions, the divergence- free con- 
dition of a current density is automatically satisfied 
through a similar logic as the CT scheme; 
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The representations for the numerical fluxes of the re- 
sistive terms in the induction equations are straightfor- 
ward, since a current density is defined at the same grid 
position as a numerical flux: 
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The numerical fluxes of the resistive terms in the energy 
equation are given as 5 
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In order to close the equation system (1) — (6), we fur- 
ther need to know a gravitational potential and the rela- 
tion between pressure and other thermodynamic quanti- 
ties. The former is done by solving the Poisson equation; 
A^> = AitGp. In Yamazakura, this is numerically solved 
by Modified Incomplete Cholesky decomposition Conju- 
gate Gradient (MICCG) method (Gustafsson 1983). For 
the latter, we adopt a tabulated nuclear equation of state 
produced by Shen et al. (1998a, b), which is commonly 
used in recent core-collapse simulations; e.g. Sumiyoshi 
et al. (2005); Murphy & Burrows (2008); Marek et al. 
(2009); Iwakami et al. (2009). To derive a pressure from 
the EOS table, three thermodynamic quantities should 
be specified: in our case, density, specific internal en- 
ergy, and electron fraction are chosen. We do not solve 
the evolution of an electron fraction. Instead they are 
assumed a function of density according to Liebendorfer 
(2005) 6 . A neutrino transport that may be important 
in the dynamics of a core-collapse is not dealt in the 
present simulations. Since a neutrino cooling is not con- 
sidered, only a photo-dissociation of heavy nuclei takes 
energy away from shocked matters. Nevertheless, our 

5 Numerical fluxes given here are a little different from those 
found in the original KT scheme. In the original scheme the average 
of cell center values, Bf and Bf +1 , are used while we employ the 

average of left and right-interpolated values, B^~ 1 ^ 2 and B*+i/ 2 m 

6 In Liebendorfer (2005), prescriptions not only for an electron 
fraction, but also for a neutrino stress and entropy change are sug- 
gested. In the present simulations, we only adopt the prescription 
for electron fraction. 
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core-collapse simulation without magnetic field and ro- 
tation still indicate that no explosion occurs (see § 4). 
Although a neutrino heating is also omitted, it may not 
be very important since the present computations are 
run until at most ~ 100 ms after bounce, a several factor 
shorter than the heating time scale. 

Although we adopt the Liebendorfer's prescription for 
electron fraction through a whole evolution, it is only 
valid until bounce. For example, that prescription does 
not properly reproduce a decrease in electron fraction 
around the neutrino sphere due to the neutrino burst. 
A numerical simulation done by Sumiyoshi et al. (2007), 
which deals with sophisticated neutrino physics, shows 
that a electron fraction around the neutrino sphere at 
100 ms after bounce is ~ 0.1, while it is ~ 0.3 in our sim- 
ulations. We have tested in the simulation without mag- 
netic field and rotation at 100 ms after bounce, how much 
a pressure around the neutrino sphere varies when the 
electron fraction is replaced from the the Liebendorfer's 
value to 0.1, and found that the difference is at most only 
20 %. Note that in the present simulations, an electron 
fraction may influence dynamics only through a pressure 
and sound speed, where the latter is just related to the 
strength of a numerical diffusion. 

Yamazakura has passed several numerical test prob- 
lems, which are demonstrated in Appendix. 

3. COMPUTATIONAL SETUPS 

We follow the collapse of the central 4000 km core of 
a 15 Mq star provided by S. E. Woosley (1995, private 
communication). To construct the initial condition of 
the core, the density and temperature distributions are 
taken from the stellar data. Note that the initial profile of 
electron fraction is determined using the prescription by 
Liebendorfer (2005) as mentioned above. In order that 
the collapse proceeds in the presence of strong magnetic 
field and rapid rotation, a temperature of the initial core 
is reduced as 

T(r) = T org (r) (l - -yll-j) (21) 

where r is the distance from the center of the core, T org (r) 
is an original temperature of the core, and rx is taken 
1000 km 7 . The initial internal energy distribution is ob- 
tained by the EOS table, using density, electron fraction, 
and temperature as the three parameters. Radial veloc- 
ities are initially assumed to be zero. Magnetic field and 
rotation are initially put by hand into the core. We as- 
sume that the initial magnetic field is purely dipole-like, 
and the core is either rotating or non-rotating. 

The dipole-like magnetic field is produced by putting 
a toroidal electric current of a 2D-Gaussian-like distribu- 
tion centered at (m,z) = (tuo,0), 

U{ia , z )=^/"(-W) (22) 



where f = y ' (w — wo) 2 + z 2 . The last factor is multi- 
plied to impose = along the pole. A width a is a 

7 In what follows we denote a spatial point in the polar and 
cylindrical coordinate by (r, </>) and (vu,cf),z), respectively. 




a [cm] 

Figure 2. Initial magnetic field configuration (vector) and distri- 
bution of specific magnetic energy in logarithmic scale for a repre- 
sentative model (Bs-fl-rj-oo; see text). Note that all models share 
the same initial field configuration. 

function of = arccos(z/f), defined by 

a{9) = " dec (23) 
V 1 — e 2 cos 

which traces a prolate ellipse centered at (w, z) = (mo, 0) 
with an eccentricity e and major radius r^ ec . Parame- 
ters in Eqs. (22) and (23) are put as zdq = 1000 km, 
^dec = 710 km and e = 0.5 in every computation. A pa- 
rameter jo, which determine the field strength, is given 
later. From the above electric current distribution, the 
vector potential is calculated by 

A^(m,z) (24) 

i r core f 2n r core u(^z c )co^ c 

= - / / — j -w c aw c d(j) c az c 

cJo Jo Jo R(m,z,m c ,<p c ,z c ) 

where R(m, z, tu c , C , z c ) is the distance between (m, 0, z) 
and (m c ,(j) c ,z c ), and m core = z core = 4000 km. The 
magnetic field is obtained via B = V x A. Note that, 
evaluating at a cell corner, the initial magnetic field 
automatically satisfies the divergence free condition as 
the same way in the CT scheme. The initial magnetic 
field configuration and the distribution of magnetic en- 
ergy per unit mass are shown in Fig 2. 

In each rotating model, an initial angular velocity is 
given by 

n W= "03X3- (25) 

where ro = 1000 km and Qo = 3.9 rad s -1 . 

Employing the above magnetic field and rotation, we 
study three different cases, namely, strong magnetic field 
and rapid rotation (model-series Bs-fi), moderate mag- 
netic field and rapid rotation (model-series Bm-fi), and 
very strong magnetic field and no rotation (model-series 
Bss-cJ). Parameters for each model-series are given in 
Table 1. 

In order to study effects of resistivity on the dynam- 
ics, both an ideal model and resistive models are run in 
each model-series. For the resistive models, we exam- 
ine two different values of resistivity, say 77 = 10 13 and 
10 14 cm 2 s _1 , reminding the discussion in § 1. Resistivity 
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Table 1 

Parameters and some results for the computed models. 



Model name 


E m /W a [%} 


T/W b [%} 


jo [cgs- Gauss] 


Bo,i c [G] 


Bs-Rot 


0.5 


0.5 


5.4 x 10 14 


9.7 x lO 1 ^ 


Bm-Rot 


0.05 


0.5 


1.7 x 10 14 


3.1 x 10 12 


Bss-Nonrot 


5.0 


0.0 


1.7 x 10 15 


3.1 x 10 13 



a Initial ratio of magnetic energy to gravitational energy. 
b Initial ratio of rotational energy to gravitational energy. 
c Initial magnetic field strength at the center. 



is uniform in space and time except that it is set zero in- 
side the radius of 10 km to save the computational time. 
There a magnetic field and thus a resistivity are expected 
to be unimportant due to high density. For a descriptive 
convenience, we use abbreviations 77-00, ?7i3, and 7714, re- 
spectively, for models with 77 = 0, 10 13 and 10 14 cm 2 s _1 , 
attaching after name of a model-series introduced above. 
For example, a model with strong magnetic field, rapid 
rotation, and rj = 10 14 cm 2 s _1 is referred to as model 

Bs-Q-TJi^. 

Each computation is done in cylindrical coordinate. 
Assuming the axisymmetry and equatorial symmetry, we 
take the numerical domain as {w, z) G [0 km, 4000 km] x 
[0 km, 4000 km]. Until the central density reaches 
10 12 g cm -3 , the number of numerical cells is x N z = 
320 x 320. After that the number of cells is changed to 
x N z = 720 x 720. There, the spatial width of a cell 
increases outward in each direction with a constant ratio 
of 1.0051 and 1.0056, before and after the re-griding, re- 
spectively. Both the width Aw of the inner-most cells of 
^-coordinate and the width Az of the inner-most cells 
of z-coordinate are 5 km and 400 m, before and after the 
re-griding, respectively. When the numerical cells are re- 
distributed, physical variables are linearly interpolated 
from the coarse into the fine cells. In this procedure, the 
divergence-free constraint of magnetic field is usually vi- 
olated, which stems from the poloidal components. To 
avoid this, we calculate from Eq. (24) using the distri- 
bution of jfj, right after the cell redistribution, and then 
obtain a divergence-free poloidal magnetic field. 

Each simulation is run until the shock front reaches a 
radius of 2100, 3000, and 2300 km in model series Bs-£7, 
Bra-f2, Bss-U, respectively. 

4. RESULTS 

In this section, we will present results of our compu- 
tations for each model-series separately, seeing how a re- 
sistivity affect the dynamics of a magnetized supernova. 
Particular attentions are paid to the explosion energy, 
magnetic field amplification, and the aspect ratio of the 
eject a. 

Before proceeding to the main results, we here describe 
the dynamical evolution in the simulation without mag- 
netic field and rotation. Soon after the start of computa- 
tion, the core has a negative radial velocity everywhere, 
and collapses towards the center. A bounce occurs at 
t = 133 ms due to nuclear force, and a shock wave is 
generated. The shock wave propagating outward first 
stalls around r ~ 200 km, but it starts gradual expan- 
sion around 165 ms. Afterward, the shock surface alter- 
natly expands and shrinks. We followed the evolution 
until t = 350 ms (217 ms after bounce) during which 
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Figure 3. Evolutions of the total (black line), internal (magenta), 
gravitational (gray), positive kinetic (orange), and negative kinetic 
(cyan) energy integrated over the whole numerical domain for the 
simulation without magnetic field and rotation. The total and 
gravitational energy are multiplied by —1. 

the maximum shock position is ~ 800 km. In this way, 
our model without magnetic field and rotation does not 
result in a stalled shock as seen recent core-collapse sim- 
ulations (see e.g. Fig. 2 of Nordhaus et al. (2010)), which 
may be because we only consider a photo-dissociation of 
heavy nuclei but no neutrino cooling as cooling processes. 
Fig. 3 shows the evolutions of total, internal, gravita- 
tional, positive kinetic, and negative kinetic energies in 
the simulation without magnetic field 8 . This indicates 
that a part of fluid has a positive radial velocity. Non- 
theless, we found that an estimated explosion energy is 
very small; less than ~ 10 48 erg and sometimes exactly 
zero. We assume that a fluid element is exploding if the 
total fluid energy plus the gravitational potential energy 
at its position, and the radial velocity are both positive, 
i.e. e + pv 2 /2 + B 2 /Stt + p$ > and v r > 0. Then 
the explosion energy is obtained by the sum of the total 
energy of fluid elements that fulfill the criterion plus the 
gravitational potential energy for the exploding fluid. We 
calculate the latter by ^ e x P ,grv = J v [p^exp/2 + p$]dV : 
where <£ e x P is the gravitational potential due to the ex- 
ploding fluid and 4> is that due to the other fluid, and 
$ = $exp + Gravitational potentials $ e xp and are 
obtained by solving a Poisson equation with only the 
mass of the exploding fluid and non-exploding fluid, re- 
spectively. The fact that the explosion enrgy is quite 
small as mentioned above implies that most or all fluid 
elements including those with a positive radial velocity 
do not fulfill the criterion. It seems reasonable to assume 
that the model without magnetic field and rotation does 
not explode. 

The error in the total energy conservation of the system 
is 51 % at the end of the simulation. We found that 
the error in the total energy conservation is 24-33 % at 
the end of the all simulations involving magnetic field, 
except that it is 14 % in a different resolution run for 
model Bs-Q-r]i4 described in § 4.1.5. In § 4.1.5, we will 
discuss whether these errors are problematic for results 
presented in this paper. 

4.1. Strong Magnetic Field and Rapid Rotation — 
Model Series Bs-Q 



8 Here, the terms "positive kinetic energy" and "negative kinetic 
energy" mean the kinetic energy associated with fluid elements 
with a positive and negative radial velocity, respectively. 
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Figure 4. Evolutions of the total (black lines), internal (magenta), gravitational (gray), rotational (red), positive kinetic (orange), negative 
kinetic (cyan), poloidal magnetic (blue), and toroidal magnetic (green) energy integrated over the whole numerical domain. The total and 
gravitational energy are multiplied by —1. The solid and dashed lines are drawn for model Bs-Q-n-oo and Bs-f2-?7i4, respectively. The 
dotted lines, which are almost identical with the dashed lines, are for a different resolution run for model Bs-fl-ni4 (see § 4.1.5) 



We start from briefly describing the dynamical evolu- 
tion in model Bs-Cl-rj-oo. In this model, a rotation ham- 
pers the collapse, and a bounce occurs at t = 143 ms, 
10 ms later than in the case without magnetic field and 
rotation. During the collapse, the core is largely spined 
up accompanied with an increase in the degree of dif- 
ferential rotation: Just after bounce, a rotational period 
reaches ~ 1 ms in a considerable part inside the radius 
of 50 km, while it is initially at least ~ 1 s. A magnetic 
field, which initially plays little role, is greatly amplified 
by compression during collapse. In addition, the differen- 
tial rotation winds poloidal magnetic field-lines, and the 
toroidal component of magnetic field is largely generated 
around the time of bounce. An outward matter motion 
driven by bounce first decelerates, losing the energy due 
to a photo-dissociation of heavy nuclei, but accelerates 
again helped by a magnetic pressure of toroidal field and 
centrifugal force. As a result, a strong eruption of mat- 
ter occurs preferentially along the pole. The above dy- 
namical sequence can be followed in Fig. 4 in terms of 
energetics (see thick lines): i.e. a decrease of the grav- 
itational energy results in an increase of the rotational 
and magnetic energy; the rotational energy is partially 
converted into the toroidal magnetic energy; then the 
toroidal magnetic energy is consumed to boost the pos- 
itive kinetic energy. The top panels of Fig. 5 shows the 
distributions of velocity and magnetic field at 164 ms 
(21 ms after bounce) for model Bs-fi- 77-00. A fast mass 
eruption (v r > 5 x 10 9 cm s _1 ) is seen notably around 
the pole, where the ratio of a magnetic pressure to matter 
pressure is large. This also implies that magnetic force 
plays an essential role for a fast mass eruption. Note 
that the dynamical features described here is quite simi- 
lar to those found in previous works that employ similar 



strengths of magnetic field and rotation (see e.g. Yamada 
& Sawai (2004), Takiwaki et al. (2009)). 

In a resistive model B s-fJ-r/u, the evolution proceeds 
in qualitatively similar way to the ideal model Bs-f}- 
rj-oc. However, outgoing velocities are relatively slow 
compared with the ideal model, which result in a smaller 
shock radius at a same physical time (compare the left 
panels of Fig. 5). The right panels of Fig. 5 implies that 
this is due to a less strong magnetic pressure in model 
7714. As easily expected, a magnetic field amplification by 
differential rotation is ineffective under the presence of re- 
sistivity. This means that the rotational energy cannot 
be spent efficiently as an energy source for the explosion. 
Indeed, it is observed in Fig. 4 that the rotational en- 
ergy in model 7714 decreases slowly compared with that 
of model 77-00 as well as both the positive kinetic and 
magnetic energy increase slowly. 

4.1.1. Explosion Energy 

Below, we will see the effect of resistivity on the explo- 
sion energy together with the detailed mechanism of the 
explosion. Fig. 6 shows the evolutions of the explosion 
energies, E exp , in model-series Bs-Sl. It is found that a 
larger resistivity results in a smaller explosion energy. 

As a preparation for detailed analyses, we consider di- 
viding the volume inside the shock surface into the two 
parts, say, the eruption-region and the inf all-region. The 
definition of these parts are as follows. First, the volume 
inside the shock surface is equally cut up with respect 
to into 30 volume segments with AO = 3° opening an- 
gle. The eruption-region is defined by the sum of the 
segments whose integrated radial momentum is positive, 
whereas the infall-region is by the sum of those having 
the negative radial momentum. For example, in each left 
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Figure 5. Left panels: Distributions of radial velocity magnitude (color) and velocity direction (vectors). Right panels: Distributions of 
ratio of magnetic pressure to matter pressure, Pm/p, in logalithmic scale (color), and magnetic field direction (vectors). These figures are 
depicted at t = 164 ms for model Bs-fl-n-oo (upper panels) and Bs-fl-ni4 (lower panels). 
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Figure 6. Evolutions of the explosion energies in model-series Bs- 
Q. The red, green, and blue-solid lines are for model ?7— oo, ??13, and 
7714, respectively. The blue-dotted line is for a different resolution 
run for model Bs-f2-?7i4 (see § 4.1.5). For the definition of an 
explosion energy, see text. 

panel of Fig. 5, the infall-region appears in the vicinity of 
the equator (6 > 70°) where a negative momentum dom- 
inates over a positive one, while the other part inside the 
shock surface corresponds to the eruption-region. 

To know what causes the differences in the explosion 
energy, it may be helpful to compare between model 



B s-Q-rj- oc and B s-Q-rji4 by the individual acceleration 
terms in the r-component of the equation of motion writ- 
ten on the frame rotating around the pole with the an- 
gular velocity 



D'v r 
Dt 



1 dp d<S> 
p dr dr 
+ft 2 rsin 2 <9, 



— UeB 4 
pc 



U B e) 



(26) 



where D' / Dt denotes the Lagrangian derivative on the 
rotating frame. In the r.h.s. of Eq. (26), each term rep- 
resents, from left to right, the acceleration due to a pres- 
sure, gravity, magnetic field, and rotation. In compar- 
ing the accelerations, we take an angular average in the 
eruption-region and a time average during t = t x -t 2 ms, 
which are defined by 



(a)(r) 



f ap sin 6d6 

Jerup ' 



f p sin OdO 



dt / [t 2 - h] , (27) 



where each term in the r.h.s. of Eq. (26) is to be as- 
signed to a. Fig. 7 shows the radial distributions of accel- 
erations, angularly averaged in the eruption-region and 
time averaged during t = 147-155 ms, in model r]-^ and 
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Figure 7. Radial distributions of the radial accelerations, (a), angularly averaged in the eruption-region and time averaged during 
t = 147 — 155 ms, in model Bs-f2-77_oo and Bs-f2-77i4. In each panel, the black, red, blue, and magenta lines represent the total, pressure, 
gravitational, magnetic, and centrifugal accelerations, respectively. Except for the right-top panel, the solid lines are for model 77-00, while 
the dashed lines are for model 7714. The right-top panel shows the differences in the accelerations, A (a) = (a) ?7 _ 00 — (a)r] 14 , subtracting 
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7714. The averages are taken inside ~140 km, the smallest 
shock radius at 147 ms in model 77-00 and 7714. It is found 
that, in the both models, an acceleration is almost every- 
where positive for r > 50 km (left-top panel). As shown 
in Fig. 8, a matter ejection also occurs roughly in this 
radial range, which implies that the amplitude of an ac- 
celeration is a good measure for the resulting magnitude 
of the explosion energy. There, a pressure acceleration 
alone is always smaller than a gravitational deceleration 
(left-bottom panel). The right-bottom panel shows that 
it is a magnetic and centrifugal acceleration that makes 
a matter ejection possible. 

In the right-top panel of Fig. 7, the differences in ac- 
celerations between the two models are plotted, which 



Figure 9. Distributions of specific angular moment urns in spher- 
ical mass coordinate at t = 155 ms in model Bs-f2-77_ 00 and Bs- 
£7-7714. The Black-dashed line shows the initial distribution shared 
by the two models. The magenta and cyan vertical lines show the 
mass that corresponds to r = 25 km, respectively, in model 77—00 
and 7714. 

shows that a total acceleration is averagely larger in 
model 77-00 for r > 50 km. It is likely that this causes 
a larger dE exp /dr in model 77-00 as observed in Fig. 8. 
The right-top panel also indicates that a larger total ac- 
celeration in model model 77-00 is primary due to that 
in a magnetic and centrifugal acceleration. Although, a 
pressure acceleration is smaller in model model 77-00, this 
is more than compensated by them. Therefore, we con- 
clude that a resistivity makes the explosion less energetic 
due to a small magnetic and centrifugal acceleration. 

While a smaller magnetic acceleration in model 7714 
will be simply due to a weaker magnetic field as a result 
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Figure 10. Left: Evolution of the average density for v ^> 50 km in the eruption-region in model Bs-Q-Tj— oo 

and Bs-f2-?7i4. Right: 

Evolution of the average specific internal energy for r > 50 km in the eruption-region in model 77-00 (red line), model 7714 (blue), and the 
difference between them (black). The magenta line shows the time integrated Joule heating produced in the above region for model 7714. 
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Figure 11. Evolutions of the explosion energies in model-series 
Bs-r2 in terms of the remaining rotational energy. The evolution 
proceeds roughly in the counterclockwise direction. 

of a magnetic diffusion, that of centrifugal acceleration 
is related to a less efficient angular momentum transfer 
owing to a weaker magnetic stress. In Fig. 9, the dis- 
tributions of specific angular momentums at 155 ms is 
plotted in mass coordinate. It is observed that a specific 
angular momentum in model 7714 is larger than that of 
model 77-00 inside the radius of ~ 25 km, whereas it is 
smaller outside there. This is the consequence of a less 
efficient outward angular momentum transfer in model 
7714. Then, at a large radius, a centrifugal acceleration in 
model 7714 is smaller than that in model 77-00. 

We also examined reasons for a larger pressure accel- 
eration in model 7714 observed in the right-top panel of 
Fig. 7. In the present situation, a pressure is a increas- 
ing function of density and specific internal energy. Al- 
though a pressure also depends on electron fraction, this 
dependence is very weaker than that on the above two 
quantities. Hence, it is expected that the density or spe- 
cific internal energy volume- averaged over the eruption- 
region are larger in model Ds-^-7714. By calculating these 
two average values, we found that both the average den- 
sity and specific internal energy are larger for model 7714 
(see Fig. 10). A higher density in model 7714 implies that 
the matter expansion rate is smaller, which is consistent 
with what is observed in the left panels of Fig. 5. A larger 
specific internal energy in model 7714 also may come from 
the smaller expansion rate, but may be caused by Joule 
heating too. To estimate an amount of thermal energy 
produced by Joule heating in model Ds-^-7714, the total 
heating rate of the eruption-region, / erup 47T77J 2 / (pc 2 )dm, 
is time-integrated from 147 ms when the infall-region be- 
gins to appear clearly. Then it is divided by a mass of the 



expulsion region at each instant of time, and compared 
with the difference in average specific internal energy be- 
tween the two models. The result is shown in the right- 
panel of Fig. 10, which indicates that the contribution 
form the Joule heating to the difference in the specific 
internal energy is quite small around 150 ms. Thus, it 
is not likely that a larger pressure acceleration in model 
7714 observed in the right-top panel of Fig. 7 is caused 
by Joule heating. A larger pressure acceleration seems 
to be just due to a smaller expansion rate in model 7714. 
The panel implies that only in a later phase, a larger 
specific internal energy in model 7714 may be somewhat 
contributed by the Joule heating. Note, however, that 
the thermal energy estimation made here may be crude, 
since a part of thermal energy produced in one region 
may migrate into another. 

As we have seen, an inefficiency both in a magnetic 
field amplification and an angular momentum transfer 
leads to a weaker explosion in a resistive model. From 
the energetical point of view, this corresponds to an in- 
efficiency in consuming the rotational energy as a fuel. 
Then one may think that the explosion energy of the 
ideal model and a resistive model would be similar com- 
pared in terms of the consumed rotational energy. If 
this is the case, the final explosion energy in the three 
models, after all the available rotational energy has been 
drained, would be comparable to each other. According 
to Fig. 11 this is not ture, however. In this figure the evo- 
lutions of the explosion energies in the three models are 
plotted against the remaining rotational energy. Since 
the maximum rotational energy, which is reached at the 
time of bounce, is almost same among the three models, 
E rot « 7.9 x 10 51 erg, each model will consume roughly 
the same amount of rotational energy with a same posi- 
tion in the abscissa. It is found that the explosion en- 
ergy in a resistive model is smaller than that of the ideal 
model, even though a same amount of rotational energy 
is expended. This implies that a part of the rotational 
energy is wasted in locations where the criterion for the 
explosion is not fulfilled. 

4.1.2. Magnetic Field Amplification 

In this section we analyze a magnetic field amplifica- 
tion. In Fig. 12, the angular distributions of the mag- 
netic energies per unit mass, averaged over 50 km< r < 
0.9xr s h at t = 145 ms (2 ms after bounce) and t = 160 ms 
(17 ms after bounce), are shown for model B s-Q-rj- 00 
and B 5-^-7714. The left panel indicates that the total 
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magnetic energy is relatively stronger around the pole at 
145 ms in each model, reflecting the initial magnetic field 
configuration (see Fig. 2). Turning to the right panel, it 
is found that, in each model, the contrast between the 
total magnetic energy around the pole and that around 
~ 40°-70° becomes stronger at 160 ms than at 145 ms. 
In model 7714, the strong contrast in the total magnetic 
energy at 160 ms is mainly due to that in the toroidal 
magnetic energy, while in model 77 _ oo, the contrast both 
in the toroidal and poloidal magnetic energies are respon- 
sible for that. 

To understand how the contrast is strengthened, we 
follow the evolution of magnetic energy per unit mass in 
two representative volumes V25.5 and V58.5, where Vq s is 
defined by 50 km< r < 0.9 x r sh and S - 1°.5 < < 
9 S + 1°.5. Fig. 13 shows the evolution of the average 
magnetic energy per unit mass in the two volumes. It is 
observed both in model rj-oc and 7714, that the toroidal 
magnetic energy in the both volumes increases around 
150 ms, and keeps an almost constant value afterward. 
The increase rate is higher in volume V25.5. That is the 
contrast in the toroidal magnetic energy is strengthened 
around 150 ms, and is kept afterward. What is also found 
is that in model 77-00, the poloidal magnetic energy in 
volume V25.5 increases during t ~ 155-160 ms, while that 
in volume V58.5 does not vary very much. It seems that 
this makes the strong contrast in the poloidal magnetic 
energy in model 77-00 at 160 ms shown in the right panel 
of Fig. 12. 

In order to study the amplification mechanisms of mag- 
netic field, we write down the evolution equations of the 
average magnetic energies per unit mass in volume Vq 3 
with mass M, <? {r> M} = [/ (B\ r ^ M /^)dV]/M: 



{r,6,cf>} 



^ " — ^{r,e, 0},adv + ^{r,6, </>},cmp + ^{r,6>,0},shr 

,</>}, rst + ^{r,6>,0},VM, (28) 

where the terms in r.h.s. mean, from the left, the change 
of S{ r ^e,<j)} due to an advection, compression, velocity 
shear along magnetic field lines, resistivity, and the vari- 
ation in the volume and mass. They are defined by 
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where v sr f is a surface velocity of volume Vq s 9 . Note that 
the resistive terms are written for a constant resistivity. 
By evaluating these terms we can see which factors are 
essential for the magnetic field amplification in a given 
volume. The results are shown in Fig. 14 and 15. 

We first focus on the amplifications of the toroidal 
magnetic energies in model 77-00- In volume V25.5, the 
toroidal magnetic energy is primarily amplified by an ad- 
vection (right-top panel of Fig. 14). We found that an 
radial advection dominates over an angular advection, 
i.e. the amplification is due to an outward advection of 
a large toroidal energy at small radii. As expected, a 
velocity shear along poloidal magnetic field-lines, viz. a 
winding due to a differential rotation, also substantially 
contributes to the amplification. In volume F58.5, the 
toroidal magnetic energy is also amplified due to an ad- 
vection and shear (see the right-bottom panel of Fig. 14), 
but the amplitude of each term is much smaller than in 
volume V25.5, which seems due to a priori weaker mag- 
netic field. 

The amplification of the radial magnetic energy in vol- 
ume 5 is also dominated by a radial advection. The 
shear term, which is far smaller than that of the ad- 
vection term, seems to also play an important role after 
~ 160 ms, since it is comparable to a total $ (see the left- 
top panel of Fig. 14). As seen in the left-bottom panel of 
Fig. 14, the amplification mechanism of radial magnetic 
energy in the volume V58.5 is rather complex, contributed 
by several terms. As in the case of a toroidal magnetic 
energy, an amplification of a radial magnetic energy is 
also weaker in volume than in volume ^25.5- 

Since the present model-series involves a rotation, the 
magnetorotational instability (MRI) may occur and may 
play an important role in a magnetic field amplification 
(Balbus & Hawley 1991; Akiyama et al. 2003; Masada 
et al. 2006). Signs for MRI growth is found in some of 
past MHD core-collapse simulations initially assuming a 
magnetar-class magnetic field (Yamada & Sawai 2004; 
Takiwaki et al. 2004; Obergaulinger et al. 2006; Shibata 
et al. 2006). Also, Obergaulinger et al. (2009) carried out 

9 These definitions are different from commonly used ones, where 
the terms on the r.h.s. of the induction equations, dB/dt = 
— (v • V)B + (B • V)v — B(V • v), are interpreted as the changes 
of magnetic field due to advection, shear, and compression (e.g. 
Brandenburg &; Subramanian 2005). We carry out possible cancel- 
lations between these terms. 
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Figure 12. Angular distributions of magnetic energies per unit mass averaged over 50 km< r < 0.9 x r s h at t = 145 ms (2 ms after 
bounce; left) and t = 160 ms (17 ms after bounce, right) in model Bs-fl-rj-oo (solid lines) and Bs-fl-rji^ (dashed lines). The dotted lines 
are for a different resolution run for model Bs-^-?7i4 (see § 4.1.5). 
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Figure 13. Evolutions of average magnetic energies per unit mass of volume V25.5 (left) and V67.5 (right). The graphs are plotted from 
t = 145 ms, 2 ms after bounce. The solid lines are for model Bs-f2-77_oo while the dashed lines are for Bs-Q-r]i4. The dotted lines are for 
a different resolution run for model Bs-^-7714 (see § 4.1.5). 



simulations of MRI with a rather weak initial magnetic 
field, which is still stronger than that of ordinary pul- 
sars by an order of magnitude, using a local simulation 
box, and found that the MRI exponentially amplifies the 
seed magnetic field. Note however that the effect of ac- 
cretion is not considered in their local box simulations. 
According to Foglizzo et al. (2006), the neutrino-driven 
convection in the gain region can be stabilized or slowed 
down by accretion. This may also hold for the MRI in 
the post-shock region. 

We investigated whether, in model r]-^ and r/14, there 
emerges a region where the criterion for the MRI (Bal- 
bus 1995) is satisfied and the growth timescale is short 
enough. In the left-top panel of Fig. 16, we plot the dis- 
tribution of a MRI linear growth timescale, roughly esti- 
mated by 47r/\u7dQ/dw\, in 0-r plane at 160 ms (17 ms 
after bounce) for model 77-00- It is shown that in a con- 
siderable part for < 40°, the growth timescale of MRI 
is a few ms to 10 ms, while in a larger the growth 
timescale is averagely longer. We found that, in the 
above part, the growth timescale of ~10 ms is kept af- 
ter bounce (t = 143 ms) until the end of the simulation. 
However, the right-top panel of Fig. 16 does not shows 
MRI-like field-line bending as observed in some models 
in Yamada & Sawai (2004) (model MF3 and MF8). This 
may be because the present model leads to a stronger 
matter eruption in the radial direction than those mod- 
els of Yamada & Sawai (2004) because of a rather mild 
initial rotation speed 10 . In the present model, field-line 



bending produced by the MRI may become invisible due 
to a dominant radial flow, and the absence of that will 
not necessarily mean the non-operation of the MRI. Note 
that the absence of field-line bending seems not due to a 
poor spacial resolution: A field-line bending is observed 
even in the computations of model-series Bm-1], in which 
a magnetic field is averagely weaker than in the present 
case and thus the resolution for capturing MRI is poorer. 
In the present case, the fastest growing wave length is re- 
solved everywhere with several 10 numerical cells (see the 
bottom panels of Fig. 16). Although the required number 
of numerical cells to capture one wave length depends on 
numerical scheme, several 10 cells will be sufficient. 

The growth of the MRI will lead to an increase of the 
shear term in Eq. (28), since the MRI produce a shear 
of velocity along a magnetic field-line. The left panel 
of Fig. 14 shows that the shear term in volume V25.5 
becomes relatively large after ~ 160 ms (~ 20 ms after 
bounce). Given that the MRI growth timescale there 
is ~ 10 ms, it may be reasonable to consider that the 
increase of the shear term is due to the operation of MRI. 
Even if this is the case, however, it is unlikely that the 
MRI greatly amplifies a magnetic field, since the radial 
magnetic energy is nearly constant after ~ 160 ms (see 
left panel of Fig. 13). The MRI seems at best to keep 
the strength of magnetic field in volume V25.5. Since the 
fastest MRI growth timescale is comparable at each for 
< 40°, the situation will be more or less similar in this 
angular range. 



It is known that a mildly rapid rotation, T/|W| ~ 0.5 %, is 



favorable for an energetic explosion (see Yamada & Sawai (2004)) 
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Figure 16. Left-top panel: Distribution of MRI growth timescale in logarithmic scale in 9-r plane. White-colored regions include a location 
that is stable against the MRI. Right-top panel: Poloidal magnetic field vectors on top of a color map of poloidal magnetic field strength. 
Left-bottom panel: Distribution of fastest growing wave length, A* ~ 27tca/^, where ca is the Alfven velocity of a poloidal magnetic field. 
Right-bottom panel: Distribution of fastest growing wave length divided by numerical cell size, A*/ \J (Ax) 2 + (Az) 2 . Yellow-colored region 
includes MRI-stable locations. These figures are depicted for model Bs-fl-rj-oo at t = 160 ms. 
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Figure 17. Same as Fig. 16 but for model Bs-Q-r]i4 at t = 160 ms. 
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Figure 18. Distribution of Rmri = T res / T mri 3 m logarithmic scale in 0-r plane in model Bs-fi-7714 at t = 160 ms. 
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Meanwhile, in the angular range of 6 > 40°, even the 
retention of magnetic field strength by the MRI will be 
modest, because a growth timescale reaches ^ 10 ms 
only in some limited locations (see the left-top panel 
of Fig. 16). Although, in volume F58.5, the shear term 
helps to increase or keep the radial magnetic energy (the 
left-bottom panel of Fig. 14), this begins long before the 
MRI is expected to grow. Only in later time (e.g. after 
~ 160 ms), the shear term might contain some modest 
contribution from the MRI. 

The amplification mechanism in model 7714 is qual- 
itatively similar to that in model rj- 00, although the 
amplitude of each $ is usually smaller due to resistiv- 
ity (see Fig. 15). The left-top panel of Fig. 17 plotted 
for t = 160 ms shows that in a considerable part for 
6 < 45°, the MRI growth timescale is 10 ms to a few 
10 ms, which is kept from the time of bounce (143 ms) 
until the end of the simulation. Fig. 18 displays the dis- 
tribution of resistive timescale divided by MRI growth 
timescale, Rmri = fres/^MRi, at 160 ms. This indicates 
that the growth of the MRI is more or less hampered 
by a resistivity. Especially, in the vicinity of the pole, 
there is almost no chance for the MRI growth. We found 
that, the growth of MRI is possible through the simula- 
tion in the angular range of 20° < < 30°, although a 
growth rate will be somewhat decreased by a resistivity in 
some locations (blue-colored regions for 20° < 6 < 30° in 
Fig. 18). In this model, the fastest growing wave length 
is resolved with several 10 numerical cells everywhere ex- 
cept for r < 20 km (see bottom panels of Fig 17). 

The left-top panel of Fig. 15 also shows that the shear 
term in volume V25.5 becomes large after ~ 160 ms. By 
the same discussion done in the above, this is possibly 
due to the MRI, but not important for a magnetic field 
amplification. Again, the MRI at best may contribute to 
keep the radial magnetic energy in model 7714. 

In light of the above discussion, it can be said that 
the contrast in magnetic energy per unit mass over ob- 
served both in model 77-00 and 7714 is strengthened or 
kept by an outward advection of magnetic energy, wind- 
ing of poloidal magnetic field-lines, and possibly by the 
MRI, all of which efficiently occur in a small-# region. 

4.1.3. Aspect Ratio 

From Fig. 5 it is found that the shape of the shock 
surface is prolate both in model B s-fl-r]-^ and model 
B 5-O-7714, in which the latter shows a less prolate fea- 
ture than the former. Defining the aspect ratio by the 
maximum- z position of ejected matters, z e j,max, divided 
by the maximum- w position, tz7 e j max , we found that it 
exceeds two at the end of each simulation (see Fig. 19). 

We first focus on the aspect ratio in model 77-00. 
Fig. 19 indicates that the aspect ratio becomes larger 
than unity soon after bounce (t = 143 ms). Since a cen- 
trifugal force is relatively stronger around the equator, it 
hampers the collapse and weakens the bounce there. This 
is one reason for the prolate matter ejection. Figs. 20 
and 21 shows the radial distributions of the accelera- 
tions, time averaged during t = 155-165 ms, in volume 
^25.5 and V58.5, respectively. It is found that in volume 
V25.5, a matter is greatly accelerated in the radial range 
of 50 km< r <160 km (see solid line in the left-top panel 
of Fig. 20). Meanwhile, in volume V58.5, an acceleration 



of matter is averagely smaller (see solid line in the left- 
top panel of Fig. 21). This implies that an acceleration 
of matter is larger for a smaller 0, which causes a fur- 
ther increase of the aspect ratio. Comparing the bottom 
panels of the two figures, a larger magnetic and centrifu- 
gal acceleration in volume V25.5 seems responsible for the 
larger acceleration. Reminding that a stronger magnetic 
field leads to a larger amplitudes of these two accelera- 
tions (see § 4.1.1), it is likely that a polar ly concentrated 
distribution of magnetic energy per unit mass (§ 4.1.2) 
is essential to generate a large aspect ratio. 

Comparing the aspect ratio among the three models in 
the left panel of Fig. 19, it is the smallest in model 7714, 
while that in models 77-00 and 7713 takes similar value. 
The right panel of Fig. 19 shows that the difference comes 
from the fact that z e j,max is largely affected by resistivity, 
while tz7 e j 5max does not change so much. With this and 
the above speculation that a prolate matter ejection is 
caused by a polarly concentrated magnetic energy distri- 
bution, the smaller aspect ratio in model 7714 is readily 
understood. In Fig. 20, it is shown that a total accel- 
eration in volume V25.5 is smaller in model 7714 due to 
smaller contribution from a magnetic and centrifugal ac- 
celeration. Meanwhile, a total acceleration in volume 
V58.5 is not very different between the two models (see 
Fig. 21), because a magnetic and centrifugal acceleration 
are not as important as in volume V25.5 due to a small 
magnetic energy per unit mass. This leads to a large 
difference only in z e j,max, and thus in the aspect ratio. 

4.1.4. Diffusion and Dissipation Sites 

In this section, we will see in which sites a resistivity 
works efficiently. In Fig. 22, the distribution of magnetic 
Reynolds number R m , the ratio of the resistive timescale 
to dynamical timescale, in model Bs-O-7713 and Bs-O- 
7714 are displayed for t = 164 ms. We define the resis- 
tive and dynamical timescales, respectively, by L 2 /rj and 
\LB\/\v x B\, where the scale length of magnetic field is 
defined by L = \cB\/\4irj\. With these definitions, the 
magnetic Reynolds number is also the ratio of the size of 
the first term to second term in r.h.s. of Eq. (5). If we 
define a diffusion-dissipation site by the location where 
the magnetic Reynolds number is <100, the right panel 
of Fig. 22 shows that the diffusion-dissipation sites in 
model 7714 are, around the pole, equator, the shock sur- 
face and in the blue filaments inside the shock surface. 
In these sites, a magnetic Reynolds number is small since 
the scale length of magnetic field is short: that of toroidal 
field is short around the pole, while that of poloidal field 
is short in the other sites. Comparing the right panels 
of Fig. 5, an intense magnetic pressure dominance seen 
around the pole in model 77-00 are not found in model 
7714. This implies that a diffusion and dissipation around 
the pole are essential for producing the dynamical differ- 
ences between model 77—00 

and 7714. 

In the case of model 7713, a resistivity works efficiently 
also around the pole, equator, and the shock surface, but 
the volumes of these sites are much smaller than those 
in the former case due to a lower resistivity. Besides, 
a magnetic Reynolds number there is generally at least 
~ 30, although in some very limited regions it reaches an 
order of unity. Hence, model 7713 is rather close to the 
ideal model. With this the explosion energy and aspect 
ratio in model 7713 are not much different from those in 
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Figure 19. Evolutions of aspect ratios (left panel) and the maximum ejecta positions in z and vo (right panel) in model-series Bs-f2. The 
dotted line in the left panel is for a different resolution run for model Bs-fl-r]i4 (see § 4.1.5). 
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Figure 20. Radial distributions of the accelerations in volume V25.5, (a), time averaged during t = 155 — 165 ms, in models Bs-Q-rj-oo 
and Bs-f2-?7i4. The right-top panel shows the differences in the accelerations, A (a) = (a) ri _ 00 — (a) ?7l4 , subtracting that of model 7714 from 
that of model 77-00- Colors and styles of lines are same as in Fig. 7 
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Figure 22. Distributions of magnetic Reynolds number R m in logarithmic scale at t = 164 ms in model Bs-fl-ni3 (left) and Bs-Q-7714 
(right). 



model f]-oo- 

4.1.5. Convergences of Results 

As a final remark for models-series Bs-Sl, we mention 
the convergences of results with respect to the size of the 
numerical cells. As a representative model, we carry out 
another run for model B s-Q-7714 with a different spatial 
resolution. Keeping the total number of cells, the size 
of the inner-most cells is changed from 400 m to 200 m, 
where a cell size increases outwards both in w and z- 
direction with the constant ratio of 1.0069. With this dis- 



tribution a resolution is higher than the original one for 
z,w < 200 km, and lower otherwise. In Figs. 4, 6, 12, 13, 
and 19, the results obtained with this different resolution 
are shown for model B by dotted lines. Although 
a little deviations from the original results are observed 
in some figures, they seem not to qualitatively and even 
quantitatively affect the discussions given above. 

The error in the total energy conservation in the differ- 
ent resolution run is 14 % as mentioned at the beginning 
of § 4. That in the normal resolution run of model Bs- 
^-?7i4, 26 %, is roughly twice as large as the different 
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resolution one. Nonetheless, as we have seen, only slight 
deviations are found between results of the two resolu- 
tion run. Hence, although the total energy error of 26 % 
may be a little too large, this will not spoil the results of 
the simulation. We expect that this is also the case for 
the other models. 

4.2. Moderate Magnetic Field and Rapid Rotation - 
Model Series Bm-Q 

The dynamical evolutions in model-series is 
qualitatively similar to those found in model-series Bs-O, 
i.e. a winding of magnetic field-lines by a differential ro- 
tation increases a magnetic and centrifugal acceleration, 
which play a key role in a matter ejection. Owing to 
a weaker initial magnetic field, the explosion occurs less 
energetically than in the former model-series. The dis- 
tributions of velocity and magnetic field at t = 181 ms in 
model Bm-Q-n-oo and Bra-^-7714 are depicted in Fig. 23. 
Compared to model-series Bs-tt with a same resistivity 
(see Fig. 5), the present model-series shows a slower out- 
ward velocity and an averagely weaker magnetic pres- 
sure. As in the former model-series, a comparison be- 
tween the upper and lower panels of Fig. 23 also indi- 
cates that an matter eruption becomes weaker with the 
presence of resistivity, with which we expect a lower ex- 
plosion energy for a resistive model. 

4.2.1. Explosion Energy 

Fig. 24 plots the evolutions of the explosion energies in 
the present model-series. The resulting explosion energy 
in each model appears to be about factor three smaller 
than the corresponding model in model-series Bs-O. It 
is shown that a resistivity also decreases the explosion 
energy as in the former case. This can be understood 
same way as we have seen in § 4.1.1 for model-series Bs- 

n. 

Fig. 25 shows the radial distribution of accelerations, 
angularly averaged in the eruption-region and time aver- 
aged during 160-170 ms, according to Eq. (27). The left- 
top panel shows that an acceleration of matter mainly 
takes place around r ~ 100 km. There, although a 
pressure acceleration alone does not overcome an inward 
gravitational acceleration (see the left-bottom panel), the 
total acceleration is positive with a help of a magnetic 
and centrifugal acceleration (see the right-bottom panel). 
From the right-top panel, it is found that a total accel- 
eration there is averagely larger in model 77-00, and that 
this is primarily due a larger magnetic and centrifugal 
acceleration. In Fig. 26, the distributions of dE exp /dr at 
170 ms are shown for model 77—00 

and 7714. It is found that 
a matter eruption occurs for r > 100 km, and dE exp /dr 
is everywhere larger in model 77-00. It is expected that a 
larger total acceleration around r ~ 100 km is responsible 
for this. Thus, it is likely that again in the present model- 
series, a resistivity makes the explosion less energetic by 
decreasing a magnetic and centrifugal acceleration as in 
model-series Bs-Q. 

4.2.2. Magnetic Field Amplification 

A magnetic field amplification in the present model- 
series goes on in a similar way as in the stronger mag- 
netic field case Bs-Q. Fig. 27 shows the distributions of 
magnetic energies per unit mass, averaged over 50 km< 



r < 0.9 x r s h at t = 143 ms (2 ms after bounce) and 
t = 180 ms (39 ms after bounce), for model B777-O-77— 00 
and B 777,-^-7714. Again, it is observed that the contrast in 
a magnetic energy per unit mass over 6 becomes stronger 
at the latter time. The evolutions of the magnetic ener- 
gies per unit mass in two representative volumes V13.5 
and V55.5 are shown in Fig. 28, which indicates that the 
difference between them becomes larger from 143 ms to 
214 ms. 

Here, we also evaluated each term in the evolu- 
tion equations of magnetic energies per unit mass (see 
Eqs. (28) and (30)). The results are shown in Fig. 29 
and 30. It is found both in model 77-00 and 7714 that the 
amplification mechanisms in volume V13.5 are similar to 
those found in volume V25.5 of model-series Bs-O, i.e. the 
toroidal and radial magnetic energy are mainly amplified 
or kept due to an advection and shear (see § 4.1.2). As in 
the former model-series, we found here that the amplifi- 
cation is weaker in a volume with a lager S , viz. weaker 
in volume V55.5. 

The left-top panel of Fig. 31 shows the distributions 
of MRI growth timescale in 6-r plane for model Bm-Q- 
77-00 at t = 170 ms (29 ms after bounce). It is found 
that a MRI growth timescale is short, ~ a few ms to 
10 ms, in a considerable part for r < 20°, while for a 
larger 0, a growth timescale is averagely much longer. 
We found that a short growth timescale mentioned above 
is kept from 148 ms (7 ms after bounce) until the end 
of the simulation. In this model, the growth of MRI is 
also assured by magnetic field-line bending, which start 
appearing shortly after 150 ms. The right-top panel of 
Fig. 31, which depicts the structure of poloidal magnetic 
field for t = 170 ms, shows apparent field-line bending 
especially around the pole, where a growth timescale is 
averagely short. 

In § 4.1.2, we have seen that in both model Bs-Q-rj-oo 
and B 5-0- 7714, the shear term for volume V25.5 starts in- 
creasing around the time when the growth of the MRI 
is expected, and thus we speculated that the increase is 
caused by the MRI. In the present model, Bra-O- 77-00, 
the shear term for volume V13.5 does not behaves like 
that. As seen in the left-top panel of Fig. 29, it has a 
large negative value around ~ 155-160 ms and a large 
positive value around ~ 170-180 ms. With this rather 
complicated behavior, it is difficult to guess when the 
MRI gives a substantial contribution to the shear term. 
Nonetheless, we could at least state that the MRI does 
not play a crucial role in amplifying a magnetic field. 
This is because the radial magnetic energy has already 
grown strong enough by the period of ~ 170-180 ms, 
while only during this period of time, the shear term 
works to amplify the magnetic field (see Fig. 28 and 29). 

In model 7714, a region of short growth timescale is not 
as limited to a small 6 as in model 77-00, albeit a growth 
timescale is averagely longer, typically a few 10 ms. As 
seen in the left-top panel of Fig. 32, although a small- 
region is more favorable site for the MRI, a large-# 
region is still subjected to a fast MRI. However, as noted 
in § 4.1.2, we should take into account that the growth 
of MRI may be suppressed by resistivity. In Fig. 33, 
the distribution of Rmri, resistive timescale divided by 
MRI growth timescale, is depicted at t = 170 ms for 
model 7714. It is shown that Rmri is smaller than unity 
in the vicinity of the pole and equator, which means no 
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Figure 23. Same as Fig. 5 but for model Bm-fl-rj-oo (upper panels) and Bm-Q-r]i4 (lower panels) depicted at t = 181 ms. 
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Figure 24. Evolutions of explosion energies in model-series Bra-Q. 

growth of MRI there. We found that, for 20° < < 40°, 
Rmri is always larger than unity in a considerable parts 
of the fast MRI regions. Hence, the MRI still seems to 
work also in model 7714. In the right-top panel of Fig. 32, 
MRI- like field-line bending is also found around the pole, 
albeit less prominent than in the ideal model. 

In model 7714, the shear term for volume V13.5 be- 
comes large at late time, ~ 190 ms (the left- top panel 
of Fig. 30). Since, as well as in the former models, the 



radial magnetic energy does not largely increase after this 
time (the left panel of Fig. 28), the shear term and thus 
the MRI does not seem to very important for a amplifi- 
cation of a magnetic field, although they may play some 
role in keeping the strength of a magnetic field. 

Compared with model-series Bs-Sl, the computations 
of present model-series have poorer spatial resolution for 
capturing MRI (see the right-bottom panels of Fig. 31 
and 32). Especially, in model Bra-^-/?-^, the fastest 
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Figure 25. Same as Fig. 7 but for model Bra-f2-77_oo and Bra-Q-7714. Radial accelerations are angularly averaged in the eruption-region 
and time averaged during t = 160 — 170 ms. In the right-top panel the difference in the total acceleration is multiplied by 5. 
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Figure 26. Radial distributions of dE exp /dr at 170 ms in model Bm-^_oo and Bra-Q-7714. 
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ular distributions of magnetic energies per unit mass averaged over 50 km< r < 0.9 x r s h at t = 143 ms (2 ms after bounce; 
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growing wave length is resolved with less than 10 nu- 
merical cells in some locations even for r > 50 km, 
where a magnetic field plays a dynamically important 
role. However, since these locations are limited to small 
areas, we do not expect that the dynamics would drasti- 
cally change, if the growth of the MRI were fully resolved 
there. 

Comparing the right panels of Fig. 12 and 27, the 
distribution of magnetic energy in model Bm-^-^-oo is 
more concentrated towards small than in model Bs-fi- 
77-00. Although the two figures are depicted at different 
times, the above mentioned feature is still observed for 
the distributions compared at the same physical time 
after bounce (e.g. 37 ms after bounce for the both mod- 
els). This may be explained as follows. Because a mag- 
netic field of model Bm-0-r/_oo is weaker, a strong matter 
eruption supported by magnetic force is restricted to a 
small region in the vicinity of the pole, where a magnetic 
field is relatively strong. Then due to an amplification 
of magnetic field by a radial advection, the contrast be- 
tween a magnetic energy in the vicinity of the pole and 
the other region becomes stronger and stronger. 

4.2.3. Aspect Ratio 

The aspect ratio of a shock surface is also become 
smaller for a larger resistivity in models-series Bm-Q (see 
left panel of Fig. 34). As seen in right panel of Fig. 34, 
the difference in the aspect ratio is attained due to that 
in £ e j,max as in the case of models-series Bs-Sl. The 
same discussion as before will explain this feature, i.e. 
a stronger magnetic field and thus a larger magnetic and 
centrifugal acceleration in a small compared with in a 
large make a shock surface prolate, which is, however, 
less standout with the presence of a resistivity. A trait of 
the present model-series is a larger impact of resistivity, 
i.e. the aspect ratio is larger in model Bm-rj-oo than in 



model Bs-77-oo, while it is smaller in model Bm-rji^ than 
in model Bs-7714. 

The larger aspect ratio in model Bm-Q-rj-oo than in 
Bs-^-77-oo would be related to the distribution of the 
magnetic energy, which is more notably concentrated to- 
wards small in the former model, as discussed in § 4.2.2. 
The smaller aspect ratio in model Bm-Q-rji^ than in Bs- 
fl-r]i4 may be understood as follows. Although the dis- 
tributions of magnetic energy is similar in their shape 
between the two models (compare the right panels of 
Fig. 12 and 27), its magnitude is averagely weaker in 
model Bra-^-7714, reflecting the initial field strength, and 
a role of magnetic field in erupting matter is less impor- 
tant than in Bs-fl-r]i4. As a result, the shape of a shock 
surface is not very much affected by magnetic field, and 
then the aspect ratio becomes smaller. 

4.3. Very Strong Magnetic Field and No Rotation — 
Model Series Bss-13 

The dynamical evolutions found in model-series Bss-c3 
are qualitatively different from those in the former two 
model-series. Without rotation, neither the generation 
of a toroidal magnetic field nor an angular momentum 
transfer occurs. However, a magnetic field still seems to 
play some role. The initially strong magnetic field af- 
fects the dynamics even during the collapse. A typical 
evolution proceeds as follows. Although a total magnetic 
energy per unit mass is a priori smaller around the equa- 
tor, that of the Bq is conversely larger there. A strong 
Be attenuates a matter infall around the equator, and 
leads to a weak bounce in the lateral direction. Due to 
the weak bounce, outgoing matters are soon prevailed by 
falling matters, and then the infall-region forms around 
the equator. Meanwhile, in the other part, a bounce 
occurs strong enough to form the eruption-region, and 
with a help of a magnetic force, especially a magnetic 
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Figure 31. Left-top panel: Distribution of MRI growth time in logarithmic scale in 6-r plane. White-colored regions include a location 
that is stable against the MRI. Right-top panel: Poloidal magnetic field vectors on top of a color map of poloidal magnetic field strength 
in logarithmic scale. Left-bottom panel: Distribution of fastest growing wave length, A* ~ 27tca/^, where ca is the Alfven velocity of a 
poloidal magnetic field. Right-bottom panel: Distribution of fastest growing wave length divided by numerical cell size, A*/ \J (Ax) 2 + (Az) 2 . 
Yellow-colored region includes MRI-stable locations. These figures are depicted for model B m-Q-n-oo at t = 170 ms. 



pressure of Bq, the shock surface propagates outwards to 
reach 3000 km at the end of the simulation. Fig. 35 shows 
the distributions of velocity and magnetic field in model 
Bss-U-r]-^ and Bss-U-r]^ at t = 195 ms (40 ms after 
bounce). In each of the left panels, the formation of the 
eruption-region and infall-region is well observed. Each 
right panel shows that a magnetic pressure is mildly im- 
portant in the eruption-region. As shown in Fig. 36 the 
magnetic energies in the two models are ~ 10 51 erg. This 
implies that a magnetic field has a potential to somewhat 
boost the explosion. Although no significant difference 
between model 77-00 and 7714 can be found in Fig. 35, as 
we will see soon later, a resistivity still plays a certain 
role in the dynamics of the present model-series. 

4.3.1. Explosion Energy and Diffusion of Magnetic Field 

In Fig. 37, the evolution of explosion energies in model- 
series Bss-t3 are plotted. As shown there, the explosion 
energies are found to be ~ 10 50 erg, one order of mag- 
nitude smaller than that of a canonical supernova. The 
figure also indicates that the models involving resistiv- 
ity produce a larger explosion energy compared with the 
ideal model, which is opposite to what is found in the 



model-serieses involving rotation. 

In Fig. 38, we plot the radial distributions of dE exp /dr 
at t = 201 ms in model Bss-U-rj-oo and Bss-U-rji^, which 
shows that the explosion energy in model 7714 is grater 
than that in model 77-00 at most radial range. We exam- 
ined whether a radial acceleration in model 7714 is lager 
than in model 77-00, but found no substantial difference 
between them. We also compared the two models by the 
radial distributions of radial velocity angularly averaged 
in the eruption-region at 201 ms (see Fig. 39). It is found 
that a radial velocity is mostly larger in model 7714. This 
seems a key factor to know why the explosion energy is 
larger in model 7714. 

Then a question is why a radial velocity is larger in 
model 7714, despite comparable strengths of radial accel- 
erations. In Fig. 40, the velocity distributions in w-z 
plane at t = 201 ms are depicted for the two models 
with a little zooming towards the center. For model 
77_ 00, it is observed that matters flow from the infall- 
region into the eruption-region around r < 200 km, and 
there damps a matter ejection. In model 7714, on the 
other hand, a matter flow into the eruption region is in- 
hibited by a "positive velocity island" observed around 
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Figure 32. Same as Fig. 31 but for model Bm-fl-r]i4 at t = 170 ms. 
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Figure 33. Same as Fig. 18, but for model Bra-f2-?7i4 at t = 170 ms. 
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Figure 35. Same as Fig. 5 but for model Bss-U-rj-oo (upper panels) and Bss-lS-7714 (lower panels) depicted at t = 195 ms. 
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Figure 37. Evolutions of explosion energies in model-series Bss- 
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Figure 38. Radial distributions of dE exp /dr at t = 201 ms in 
model Bss-U-r]-oo and Bss-LS-7714. 



150 km < r < 350 km in the infall-region, and a coher- 
ent outflow of matter in the eruption-region is well kept 
even for r < 200 km. In this way, the appearance of the 
positive velocity island seems to be related to the dif- 
ference in the velocity distribution observed between the 
two models. 

We found that the positive velocity island begins to 
emerge around 170 ms (15 ms after bounce) in model 
?7i4- Fig. 41 shows the radial distributions of radial veloc- 
ities, angularly averaged in the infall-region at 175 ms, 
for the two models. It is observed in the radial range 
of 60 km< r <180 km that a infall velocity is slower 
in model 7714, which seems due to the positive velocity 
island. In Fig. 42, the radial distributions of radial ac- 
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Figure 39. Radial distributions of radial velocities angularly av- 
eraged in the eruption-region, at t = 201 ms in model Bss-U-r]- 00 
and Bss-U-r]i4. 

celerations, angularly averaged in the infall-region and 
time averaged during t = 170-175 ms, is plotted. As ex- 
pected, a radial acceleration in the infall-region is larger 
in model 7714 around a similar radial range to the above, 
90 km< r <180 km. We found that a larger pressure 
and magnetic acceleration are responsible for this ac- 
celeration superiority (see right panel of Fig. 42). We 
speculate that a larger pressure acceleration in model 
7714 will not be caused by the Joule heating, since the 
heating timescale in the above range is too long, ~1- 
10 s, to produce an extra pressure. It is more likely that 
the accumulation of falling matter stemmed by magnetic 
acceleration makes a pressure in the positive velocity is- 
land larger. It seems that a magnetic acceleration is the 
primary factor for the formation of the positive velocity 
island. 

In the absence of toroidal magnetic field, a radial mag- 
netic acceleration is written as 



Be 

a B = — 
P 



dB e 1 dB r B e 



dr 



d0 



(30) 



The radial distributions of magnetic field in the infall- 
regions at 175 ms plotted in the left panel of Fig. 43 
indicates that B r is far grater than Bq in the above radial 
range of 90-180 km, and thus the second term in r.h.s. of 
Eq. (30), a part of magnetic tension, is dominant. The 
left panel of Fig. 43 also shows that, in the above radial 
range, Bq is larger in model 7714 than in model 77-00, while 
B r is smaller in model 7714. We found that the superiority 
in Bq overwhelms the inferiority in B r , i.e. B r (r)BQ(r) 
is averagely larger in model 7714 in the above radial range. 
To put it simply, a larger Bq in in model 7714 is crucial to 
the formation of the positive velocity island. 

The radial distribution of Bq in model 7714 shown in the 
left panel of Fig. 43 compared with that of model 77-00 
invokes a magnetic field diffusion. Indeed, it is found 
that a magnetic Reynolds number is smaller than unity 
almost everywhere for r < 85 km in the infall-region 
(see the right panel of Fig. 44), viz. a diffusion effec- 
tively advects a magnetic field outwards against matter 
infall. The right panel of Fig. 43 shows that Bq in model 
77-00 is highly confined into the central region of 30 km 
radius, while the ^-distribution in model 7714 is rather 
gentle. This implies that in the latter model, a resis- 
tivity effectively diffuses Bq out of the central region. In 
model 77-00, the gradient of Bq suddenly becomes shallow 
around r = 30 km. Then with a resistivity, the incoming 
diffusion flux there, t]\7Bq, would be larger than the out- 
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Figure 40. Distributions of radial velocity magnitude (color) and velocity directions (vectors) at t - 
and Bss-U-r]i4 (right). 
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Figure 41. Radial distributions of average velocities in the infall-region at t = 175 ms in models Bss-U-rj- 



and Bss-U-rji^. 
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Figure 42. Left: Radial distributions of radial acceleration, (a), angularly averaged in the infall-region and time averaged during t 
170 — 175 ms, for Bss-U-rj-oo and Bss-U-r]i4. Right: The radial distributions of differences in the averaged radial accelerations, A (a) 
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going one, which would result in a increase of Bq outside 
the radius of 30 km. In this way, a magnetic field diffused 
out from a deep inside of the core suppress a matter in- 
fall and responsible for a larger explosion energy in model 

A similar mechanism also works in model 7713. As 
shown in the left panel of Fig. 44, diffusion of magnetic 
field in the infall-region also occurs well in this model. 
However, due to a smaller resistivity, diffusion sites are 
limited to a smaller volume compared with model 7714. 



Nonetheless, the explosion energy in model 7713 is com- 
parable to that of model 7714 at the end of the simulations 
(see Fig. 37) 

4.3.2. Aspect Ratio 

The aspect ratios in model-series Bss-c3 are also larger 
than unity, i.e. the shape of ejecta is prolate, but are 
not as large as in the rotating cases (see Fig. 45). Since 
the initial magnetic field assumed in the present model- 
series is very strong, it affects the dynamics even before 
bounce. Due to the dipole-like configuration, a magnetic 
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Figure 44. Distributions of magnetic Reynolds number in logarithmic scale at t - 
The volume below the gray line corresponds to the infall region. 

field hampers a matter infall especially well in the lat- 
eral direction. This causes a weak bounce around the 
equator, and thus a prolate matter ejection. 

The difference from the former model-series is that the 
aspect ratios do not grow drastically after bounce (see 
Fig. 45). We expect that this is caused by a roughly flat 
angular distribution of magnetic energy per unit mass in 
the eruption-region as indicated in Fig. 46 for model Bss- 
13-T]- oo and Bss-U-rji^, reminding that a polar ly concen- 
trated magnetic field distribution is essential for a large 
aspect ratio in the rotating cases. 

We found that the change from the rather polarly con- 
centrated distribution of magnetic field at the beginning 
to a roughly flat angular distribution occurs during the 
collapsing phase, i.e. before bounce. Fig. 47 shows 
the angular distributions of § r , radially averaged over 
50 < r < 1000 km at 122 ms (33 ms before bounce) in 

model Bss-13-rj-oo. It is shown that a total S r is larger 
around the equator, and a compression of B r in the 6- 
direction plays an important role. Since Bq is a priori 
strong around the equator due to the assumed magnetic 
field configuration, an infall matter there is a little de- 
flected into non-radial direction by a Maxwell stress of 
B r Bo/4:7r, viz. ^-component of velocity is generated. In 
this way, B r is compressed into the ^-direction, which 
leads to an increase of B r preferentially around the equa- 
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Figure 45. Evolutions of aspect ratios in model series Bss-£3. 

tor and rather flat angular distribution of magnetic field. 

A notable point here is that a rather fiat distribution 
of magnetic field, which causes a mild value of the aspect 
ratio, is established before bounce when a resistivity does 
not work well because of a large scale height of magnetic 
field. This is why the aspect ratio for a different resistiv- 
ity results in a similar value. 

4.4. Numerical Resistivity 

In the preceding sections, we have seen how a resis- 
tivity impacts on the dynamics of a strongly magnetized 
supernova. However, we should remember that the nu- 
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Figure 46. Angular distributions of magnetic energies per unit 
mass averaged over 50 km< r < 0.9 x r s h at t = 180 ms (25 ms 
after bounce). The solid and dashed lines are for models Bss-U- 
77_oo and Bss-U-r]i4, respectively. 
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Figure 47. Angular distributions of <f r , radially averaged over 
50 < r < 1000 km (see Eq. (30)) at 122 ms (33 ms before bounce) 
in model Bss-13-t]-oo- The black line shows total <# r , while the red- 
solid, red-dashed, green-solid, green-dashed, blue, and cyan lines 
represent the contribution to S r from the radial advection, angu- 
lar advection, radial compression, angular compression, shear, and 
change in mass, respectively. 

merical results are not only affected by a physical resis- 
tivity but may also be influenced by a numerical resistiv- 
ity. Here, we estimate a possible impact of a numerical 
resistivity on the dynamics. 

First, we directly compared the strengths of a physical 
resistivity and a numerical resistivity. We consider di- 
viding a numerical flux into two terms, the advection 
term and numerical diffusion term, where the former 
may account for an advection of a physical flux, while 
the latter for a numerical diffusion. For example, in the 
equation of numerical fluxes F x in Eqs. (10), the first 
and second term in r.h.s. are the advection term and 
numerical diffusion term, respectively. Accordingly, the 
ideal part of a numerical flux vector for the induction 
equations F, which are defined by a combinations of 6th 
- 8th components of numerical flux vectors F x an F y 
(see Eqs. (15)), are also divided into the advection term 
and numerical diffusion term. Then, the total numerical 
flux vector for the induction equations can be written as 
E = ^"adv + ^nd + G, where subscripts "adv" and "nd" 
stand for advection and numerical diffusion, respectively. 
Recall that a term G arises due to a physical resistivity. 

Now, to compare the strengths of a physical resistivity 
and a numerical resistivity we evaluate £ = |G|/|-F n d|> m 
which a physical resistivity is dominant for £ ^> 1 while 
a numerical resistivity is dominant for £ <C 1. The dis- 
tribution of ( for model Bs-ft-r]is at 164 ms is shown in 



Fig. 48. This implies that a physical resistivity dominates 
over a numerical resistivity in a considerable volume in- 
side the shock surface, but in some locations a numerical 
resistivity is not negligible. 

Next, we examined how much a numerical resistiv- 
ity affects the dynamics. This is done by evaluating 
#m,num = | ^adv I / 1 ^nd | , which may corresponds to a 
magnetic Reynolds number for a numerical resistivity. 
Even if a numerical resistivity dominates over a physi- 
cal resistivity, it does not make substantial impacts on 
the dynamics provided -R m ,num ^ 1- Fig. 49 shows the 
distribution of i? m ,num for model Bs-^-7713 at 164 ms. 
It is found that i? m ,num is much larger than unity al- 
most everywhere except around the shock surface, where 
a steep variation of a magnetic field exists. Meanwhile, 
we have seen in § 4.1.4 that a magnetic Reynolds number 
of physical resistivity in model B 5-^-7713 is ~ 30 around 
the equator (see the left panel of Fig. 22). These facts 
imply that a numerical resistivity is mostly too small to 
influences the dynamics, while a physical resistivity af- 
fects the dynamics, albeit a little. Thus, it seems that 
a physical resistivity "effectively" dominates over a nu- 
merical resistivity in this model. 

From the above discussion, it is likely that a compar- 
ison between R m and R m ,num is more meaningful than 
that between a physical resistivity and numerical resistiv- 
ity themselves in order to assess an influence of a numer- 
ical resistivity. Then, we also compare R m and i? m ,num 
for model Bra-^-7713 and Bss-13-rjis. In model Bm-^-7713, 
^m,num is much larger than unity almost everywhere ex- 
cept around the shock surface, while a relatively small 
R 1-10 is found around the equator (see Fig. 50). 
Hence, a physical resistivity also effectively dominates 
over a numerical resistivity in this model. In model Bss- 
£3- 7713, a mildly small i? m ,num ~ 10 appears not only 
around the shock surface but also other locations (see 
the left panel of Fig. 51). This means that a numerical 
resistivity somewhat affects the dynamics in this model. 
However, the right panel of Fig. 51 shows that a physical 
resistivity is more influential than a numerical resistiv- 
ity especially at small radii. Recalling that the essential 
role of a resistivity in this model is to diffuse a magnetic 
field from a small radius to a larger radius, a physical 
resistivity seems to play primary role to yield a different 
result for Bss-U-rjis from the ideal model. For confir- 
mation, we carried out the above analysis in the models 
with 77 = 10 14 cm 2 s _1 , and found that the "effective" 
dominance of a physical resistivity over a numerical re- 
sistivity is more pronounced than the counterpart models 
with rj = 10 13 



cm 2 s 



5. DISCUSSION AND CONCLUSION 

We have done MHD simulations of CCSNe to under- 
stand roles of a turbulent resistivity in the dynamics. As 
a result, we found that a resistivity possibly has a great 
impact on the explosion energy, a magnetic field ampli- 
fication, and the aspect ratio. How and how much a re- 
sistivity affects on these depend on the initial strengths 
of magnetic field and rotation together with the strength 
of a resistivity. 

In the rotating cases (model-series Bs-fl and Bm-Q), 
a resistivity makes the explosion energy small. This is 
mainly ascribed to a small magnetic and centrifugal ac- 
celeration due to an inefficiency in magnetic field am- 
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Figure 48. Distribution of £ = |G/-F n d| in logarithmic scale for 
model Bs-f2-?7i3 at 164 ms. 
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Figure 49. Distribution of i? m ,num = l^advl/l^ndl m logarith- 
mic scale for model Bs-fl-rjis at 164 ms. 

plification and angular momentum transfer, respectively. 
Meanwhile, in the case of very strong magnetic field and 
no rotation (model-series Bss-U), a resistivity enhances 
the explosion energy. In the ideal model, an inflow of 
negative radial momenta from the infall-region into the 
eruption-region inhibits a powerful mass eruption. With 
a resistivity, a magnetic field in the infall-region diffuses 
outwards from a deep inside of the core to counteract the 
negative momentum invasion. 

The ideal models involving rotation shows a polarly 
concentrated distribution of magnetic energy per unit 
mass. Since an initial magnetic energy per unit mass 
is somewhat larger around the pole than the equator, 
an amplification preferentially occurs there, which makes 
the contrast stronger. We found that the main mech- 
anisms of amplification around the pole is an outward 
advection of a magnetic energy from a small radius and 
the winding of poloidal magnetic field-lines by differen- 
tial rotation. In a resistive model, an amplification oc- 
curs qualitatively similar way to the corresponding ideal 
model, but proceeds less efficiently. As a result, the an- 
gular distribution of magnetic field becomes less concen- 
trated towards the pole. Although we found signs of a 



MRI growth in these models, it seems not to play sub- 
stantial role in amplifying a magnetic field. The MRI 
will at best help to keep the strength of a magnetic field. 
We should note here that in our 2D-axisymmetric sim- 
ulations, the MRI may not be followed correctly, since 
non-axisymmetric modes are also important in the MRI. 
A comparison between a 2D-axisymmetric and 3D sim- 
ulations in the behavior of the MRI will be made and 
presented elsewhere in the near future. 

Every model involving a rotation shows a large aspect 
ratio, > 2.5, at the end of the simulation, where a model 
with a larger resistivity results in a smaller aspect ra- 
tio. The impact of resistivity on the aspect ratio is more 
standout in the moderate magnetic field case, model- 
series Bra-Q. The aspect ratios obtained in the rotation 
models are expected to grow even larger later on, since a 
radial velocity is still faster around the pole than around 
the equator at the end of the simulations. Further long 
time computations are necessary to predict the aspect 
ratio after the shock surface breaks out the stellar sur- 
face. In a model without rotation, on the other hand, the 
aspect ratio keeps mild value, ~ 1.5, through each sim- 
ulation. No significant difference is found among models 
with a different resistivity. It is found that the aspect ra- 
tios attained in the present simulations are related to the 
angular distribution of magnetic energy per unit mass. 
That is, the more a magnetic energy is concentrated to- 
wards a small the larger the aspect ratio of ejected 
matters. 

In literature, we have found eight CCSNe in that 
both the upper and lower limit of the aspect ratio 
is measured; SN1987a, SN1993J, SN1997X, SN1998S, 
SN2002ap, SN2005bf, SN2007gr, and SN2010jl. The ob- 
served aspect ratios are at most « 3; e.g. « 2 — 3 for 
SN1987A (Papaliolios et al. 1989), « 3 for 1997X and 
1.2 - 2.5 for SN1998S (Wang et al. 2001). It seems that 
aspect ratios obtained in our rotational models are too 
large compared with these observations, reminding that 
ours will grow larger later on. However, we note that the 
above sample does not necessarily contain supernovae 
that leave a magnet ar. Applying SGR/AXP birth rate 
of ~ 0.1 — 0.2 per century (Leahy & Ouyed 2007) and 
Galactic CCSN rate of ~ 0.8 — 3.0 per century (Diehl 
et al. 2006), the rate of magnetar production among all 
CCSNe is - 3 - 25 %. 

A constant resistivity that we assume in our compu- 
tations may not be natural. A resistivity arising from 
turbulent convective motions may appear only around 
convectively unstable regions and may take a different 
value at a different position. Computations implement- 
ing such a non-constant resistivity should be done in fu- 
ture works. Effects of viscosity, which is expected to arise 
together with resistivity in a turbulence, also have to be 
investigated. Although the ideal way to study effects of 
turbulent resistivity and viscosity is a direct numerical 
simulation of turbulence itself, it is currently quite unfea- 
sible. That kind of simulation may demand a computa- 
tional facility dozens orders of magnitude more powerful 
than those of the present day. 

Although we carried out 2D-MHD simulations with 
the assumptions of axisymmetry and a purely poloidal 
magnetic field at the beginning, it is well known that a 
purely poloidal and purely toroidal magnetic fields are 
both unstable against non-axisymmetric perturbations. 
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Figure 50. Distribution of i? m ,num (left) and R m (right) in logarithmic scale for model Bra-f2-?7i3 at 182 ms. 
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Figure 51. Distribution of R 



m.num 



(left) and R m (right) in logarithmic scale for model Bss-lS-7713 at 195 ms. 



According to Braithwaite (2009), the ratio of the poloidal 
magnetic to total magnetic energy, E p /E, should be less 
than 0.8 for the stability. In each of our rotational 
model, we find that E p /E < 0.8 is always satisfied af- 
ter bounce. Although the stability criterion is not ful- 
filled before bounce, the instability will not grow sub- 
stantially, since it takes at most one Alfven timescale 
until bounce. Thus, we expect that the results obtained 
in our rotating models will not change drastically even 
if we carry out 3D-MHD simulations. In fact, 3D core- 
collapse simulations with strong magnetic field and rapid 
rotation done by Kuroda & Umeda (2010) shows qual- 
itatively similar results to those of our 2D simulations. 
However, a purely poloidal magnetic field at the begin- 
ning itself may be unnatural. Also, in our non-rotating 
models, the above stability criterion is always not satis- 
fied. The toroidal component of magnetic field should be 
added to the initial condition to deal with more realistic 
situations. It may be necessary in the future work to in- 



vestigate how such an additional component affects the 
results presented here. 

In the present work, we studied the role of a resistivity 
in the limited parameter range of magnetic field and rota- 
tion. As we found that a resistivity affects the dynamics 
differently for a different strength of magnetic field and 
rotation, it may be interesting to carry out a systematic 
study with a wide range of the parameters in the future. 
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APPENDIX 



Here we present numerical results of various HD and MHD test problems calculated with Yamazakura. Among them, 
the numerical solution of ID Riemann problem, Sedov's self-similar point source explosion, Yahil's self-similar collapse, 
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Figure 52. Numerical results (red crosses) of the lD-Riemann problem on the reference solutions (black lines) at t = 0.164 
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linear wave propagation problem, magnetic field diffusion problem, can be compared with the analytic solutions. The 
others are multi-dimensional problem, and no reference solution exits. One can find numerical results of such problems 
in literature, and may check whether our results match up with those. In each test calculation, the ideal gas EOS, 
p = (7 — l)e, is used. The distribution of numerical cells are uniform unless otherwise stated. 

Riemann Problem in ID 

An lD-Riemann problem, which is also referred as a shock tube problem or Sod problem, is carried out to check 
that 1D-HD equations are properly solved. The initial set-up is as follows: 



'p\ /1.0^ 
v = 0.0 
kpJ \1.0, 



for 



x < 0.5, v 



'0.125 N 

0.0 
v 0.1 , 



for x > 0.5 



(i) 



with the computational domain x G [0,1]. The adiabatic index is set 7 = 1.4. Numerical results with 800 cells are 
shown in Fig. 52, in which a good agreement with the reference solution is found. 

Riemann Problem in 2D 

In order to test a multi-dimensional HD computation, a Riemann problem in 2D is carried out. A square numerical 
domain (x,y) G [0, 1] x [0, 1] is divided into the four parts, and constant hydrodynamical values are set in each part. 
The initial condition is 



for x < 0.5 and y < 0.5, 
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for x < 0.5 and y > 0.5, 



for x > 0.5 and y < 0.5, 
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Figure 53. Numerical result of the 2D-Riemann problem at t = 0.25 (density color map with contours and velocity field vectors). The 
computation is done with 400 x 400 cells. 
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Figure 54. Numerical results of a point source spherical explosion for density (left-top), velocity (right-top), and pressure (left-bottom) 
distribution. The right-bottom panel is shown to see, with linear scale, a sharpness of the numerical solution for pressure. In each figure, 
the red, green, and blue crosses represent solutions at t = 0.05, 0.2, and 0.4, respectively. 



for x > 0.5 and y > 0.5, 



(2) 

which is symmetric across a diagonal line y = x. An adiabatic index of 7 = 1.4 is used here. These setups are same 
with Case 12 of Liska (2003). The calculation is done with 400 x 400 numerical cells. A numerical result is shown in 
Fig. 53, which is good agreement with that of Liska (2003). 

Point Source Spherical Explosion 

Another HD test calculation we have done is a spherical point source explosion, the solution of that known as 
the Sedov's solution. An ID spherically symmetric computation is done with polar coordinate. We use non-uniform 
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Table 2 

Errors in numerical solution at the inner-most cell. 



Coordinate 


{Ar min , Atu min , Az min } & [cm] 


D-error b [%] 


V-error c [%] 


Cylindrical 


2 x 10 4 


7.28 


22.0 


Cylindrical 


4 x 10 4 


13.1 


22.9 


Polar 


4 x 10 4 


4.18 


1.22 



a Spatial resolution of the inner-most cell 
b Error of D at the inner-most cell 
c Error of V at the inner-most cell 

numerical cells with the spatial resolution of the i-th cell being Ar^ = Ar m i n c^ _1 , where Ar m i n = 10 -4 , a = 1.0056, 
and the maximum of i is 720. This cell distribution is similar to that of the MHD core-collapse simulations presented 
in the present work, where the same cell distributions are set along the w and z directions separately. Hence, we may 
be able to speculate here how much a propagating shock wave diffuses in our MHD core-collapse simulations. We put 
a static sphere of radius r = 1.0 with the uniform density of po = 1.0. A point source explosion is initiated by injecting 
the internal energy of Eo = 3.0 at the central sphere of radius r = 0.01. Numerical results are shown in Fig 54. It is 
found that a sharpness of density peak in the numerical solution is kept during propagation, despite that the spatial 
resolution is lower for a larger radius. This implies that, in our computations, a propagating shock wave is not damped 
crucially due to a numerical diffusion. 

Self- Similar Collapse 

In order to check gravity is correctly dealt in Yamazakura, we tested a self-similar collapse, the solution of that 
known as Yahil's solution (Yahil 1983). Both an ID spherically symmetirc calculation with polar coordinate and a 2D 
axially symmetric calculation with cylindrical coordinate are done. In describing Yahil's solution, a radius, density, 
and velocity are written in dimensionless form: 

R = K- 1 ^G^- 1 ^ 2 r(t -ty- 2 

D(R)=pG(t -t) 2 (3) 

V(R)=v^ 2 G^- 1 ^ 2 (t -tr-\ 

where G, 7, ft, and to are respectively the gravitational constant, adiabatic index, polytrope coefficient, and the time 
when the central density diverge. We put 7 = 1.3, n = 1.2435 x 10 15 /2 4 / 3 cgs, and to = 0.1 s. The collapse of a sphere 
with 4000 km radius is followed. In the cylindrical coordinate calculations, two different spatial resolutions are tested, 
where the resolutions of the inner-most cells are {Ati7 min , Az m i n } = 2 x 10 4 and 4 x 10 4 cm. The number of numerical 
cells are x N z = 720 x 720 in the both calculations. In the polar coordinate calculation, only one spatial resolution, 
Ar m i n = 4 x 10 4 cm, with N r = 720 are taken. 

The distributions of density and velocity when the central density reaches 2.3 x 10 14 g cm -3 are presented in Fig. 55, 
while errors at the inner most numerical cells are shown in Table 2. Note that, in our MHD core-collapse simulations, 
a bounce occurs when the central density reaches ~ 2 x 10 14 g cm -3 . Although the spherical collapse of the 15 M 
progenitor used in the present work is not very well described by Yahil's solution 11 , the results obtained here may 
indicate how much error our collapse simulations involve. From the left panel of Fig. 55 one can see that, a deviation 
from the analytic solution is large near the center in each calculation 12 . In the cylindrical coordinate calculation with 
{Ati7 min , Az m in} = 4 x 10 4 cm, in which the grid construction is same with our standard MHD collapse simulations, 
the errors both in D and V are ~ 10 % (see the second row of Table 2). Although these errors seem too large to pursue 
a proper numerical simulation, as discussed in § 4.1.5, varying the spatial resolution in a MHD collapse simulation 
does not change results very much. This implies that the errors around the center are not fatal to the extent of the 
discussions developed in this paper. 

As seen in Fig. 55 and in Table 2, the cylindrical coordinate calculations have larger errors than the polar coordinate 
calculation does, which may be due to the structure of numerical cells. However, we found that numerically evaluated 
electric currents behave better in cylindrical coordinate than in polar coordinate, and thus adopt the former in the 
present MHD simulations. 

MHD Riemann Problem 

As a basic test for the ideal-MHD part of our code, we carried out a 1D-MHD Riemann problem that are known as 
Brio-Wu Problem (Brio & Wu 1988). The initial condition is set as follows; 
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11 We found that when the central density is in the range of ~ 10 12 — 10 13 g cm 3 in a spherical collapse of the progenitor, the 
distributions of a density and velocity is roughly described by Yahil's solution. Out of this range, a deviation from Yahil solution is large. 

12 The deviations from the analytic solution around the outer boundary are a numerical artifact due to a boundary condition. 




0.01 0.1 1 10 100 100C 0.01 0.1 1 10 100 1000 

R R 



Figure 55. Results for the calculations of the Yahil's self-similar collapse with cylindrical coordinate (left) and polar coordinate (right). 
The solid lines and dotted lines are, respectively, for the analytic solution of D and V, while the colored crosses represent the numerical 
ones. For the cylindrical coordinate calculations, the numerical solutions along the pole (vj=0) are shown. 
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Figure 56. Results of the MHD-Riemann problem at t = 0.1. The red and Green crosses represent the solutions obtained by Yamazakura 
and ZEUS-2D, respectively. The both calculations are done with 800 numerical cells. 

(4) 

where the computational domain is x £ [0,1]. The initial velocities are set to be zero. The adiabatic index here is 
2.0. The results of the calculation at t = 0.1 done with 800 numerical cells are given in Fig. 56 together with those 
obtained by a famous numerical code ZEUS-2D (Stone & Norman 1992). Although a little differences can be found, 
the two results well matches. 



38 



Sawai et al. 



Linear Wave Propagation 

In § 2, we mentioned that the KT scheme adopted in our numerical code is third order in time and second order 
in space. Here, we will see whether Yamazakura solve the MHD equations fairly with these degree of accuracy, by 
carrying out a linear wave propagation problem in ID Cartesian coordinate. A right-going slow magnetosonic wave is 
adopted as a linear wave. 

According to Gardiner & Stone (2005), we put the initial state of the conserved variables as 



Uin = u + e J Rcos(27rx), 



(5) 



where Uo is the background state, e = 10 6 is the wave amplitude, and R is the right eigenvector for the right-going 
slow magnetosonic wave in conserved variables given by 
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(6) 



The background state is set as po = 1, po = 1, vq x = vo y = vq z = 0, Bq x = 1, Bq v = V% and Bq z = 1/2, with 
which the slow magnetosonic speed is c s = 1/2. The adiabatic index is set as 7 = 5/3. The computational domain is 
x G [0,1], and the periodic boundary condition is adopted at each boundary. Each calculation is run until t = 4 when 
the wave has propagated the distance of two wave length. 

In order to check the degree of accuracy in space, we first tested two calculation serieses with uniform cell construction, 
each of which contains four calculations with different grid numbers. In the first series, series A, we fix At by controlling 
the Courant-Friedrichs-Lewy (CFL) number while varying Ax. The numbers of cells N and CFL numbers v for the 
four calculations are (TV, v) = (8, 0.49/64), (32, 0.49/16), (128, 0.49/4), and (512,0.49). In the second series, series B, 
we fix the CFL number at v = 0.49 and tested the same set of the cell numbers as before. Note that in series B, 
both Ax and At are halved by doubling the cell number. The distributions of v x at t = 4 for series B are plotted 
in Fig. 57, which shows that the numerical solution well matches the analytic one for TV > 128. The error from the 
analytic solution is evaluated by the LI error norm, 
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The left panel of Fig. 58 plots the variation of LI error norm with respect to TV for series A and B. The red plots (series 
A) show that the scheme is approximately second order in space. We found that the results from series B (blue plots) 
are almost identical to those in series A, despite that At is larger in series B for a fixed TV, viz. a fixed Ax, except 
for TV = 512. This indicate that an error due to At is dominated over by an error due to Ax for v < 0.5 (the CFL 
condition for the KT scheme), and brings difficulty in properly checking the degree of accuracy in time. Nonetheless, 
we can still state that the scheme is at least approximately second order in time, by the fact that the LI error norms 
in series B, varying both Ax and At at a same rate, result in ~ TV -2 (see the blue plots in the left panel of Fig. 58). 

Next, we evaluate the degrees of accuracy of our code for cases of non-uniform cell construction. A structure of 
numerical cells we set here is such that the (2n — l)-th cell has the size L/N x (1 — a) while the 2n-th cell has the size 
L/N x (1 + a), where L is the size of the whole numerical domain and a < 1 is the parameter to represent the degree 
of non-uniformity. We tested three cases a = 10 -2 , 10 _1 , and 5 x 10 _1 . In each case, four calculations fixing v = 0.49 
with TV = 8,32,128, and 512 are done. The results are plotted in the right panel of Fig. 58. It is shown that the 
second order in space and at least second order in time are marginally kept even for 200 % size increase and 33 % size 
decrease from the neighboring cell (a = 5 x 10 _1 ). Note that a numerical cell size in our MHD-collapse simulations 
increases outwards by 0.56 %. 



Rotor Problem 

In order to check that the multi-dimensional ideal MHD equations are well solved with our code, we perform two 
test calculations. One is a so-called rotor problem to be presented here, and the other is a point source explosion with 
magnetic field that will be shown in the next section. 

The rotor problem is suggested by Balsara & Spicer (1999), which simulates the propagations of torsional Alfven 
waves. The system initially consists of a rapidly rotating cylinder of dense fluid surrounded by lighter gas at rest, 
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Figure 57. Numerical results of the linear slow magnetosonic wave propagation for series A for the distribution of v x at t = 4. The plots 
for the calculations with N = 128 and 512 are thinned for the purpose of clear illustration. 
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Figure 58. Distributions of LI error norm with respect to the number of numerical cells for the linear slow magnetosonic wave propagation 
problem with uniform cell spacing (left) and non-uniform cell spacing (right). 

which is threaded by uniform a magnetic field parallel to the x-axis. With the computational domain of (x, y) G 
[—0.5,0.5] x [—0.5,0.5], the initial condition is 

p = l + 9/(r), 

-2f(r)y/r, v y = -2f(r)x/r, v z = 0, 
5, B y = B z = 0, 



v x 
B x 
where 



(9) 



if r<0.1, 
if 0.1 < r 
if r> 0.115, 



^(0.115 -r) if 0.1<r< 0.115, 



2U/2 



The adiabatic index is 7 = 1.4. The calculation is done with 400 x 400 numerical cells. Results of the calculation are 
displayed in Fig. 59, which show good agreement with those found in Ziegler (2004). 

Point Source 3D -Explosion with Magnetic Field 

A calculation presented here differs from the spherical explosion test described above in that fluid is magnetized. 
The initial condition is given by 



10 4 if x 2 +y 2 + z 2 < r 2 , 



1 



J x — u y 



otherwise. 
v z =0, 



5^/^4^=100, B X = B : 



y 



0. 



(10) 

The adiabatic index is chosen as 7 = 5/3. We do a calculation in cylindrical coordinate with the numerical domain 
of {vo, z) G [0,0.5] x [0,0.5], taking the symmetric axis in the magnetic field direction. The number of numerical cells 
is x N z = 96 x 96. A result of the calculation is given in Fig. 60. This also shows good agreement with that of 
Ziegler (2004) done with Cartesian coordinate. 
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Figure 59. Numerical results of the rotor problem. Left and Right panels respectively show pressure and Mach number contours at 
£ = 0.15 




Figure 60. Magnetic pressure contours at t = 2.5 x 10 3 in the calculation of a point source explosion with magnetic field. 



Magnetic Field Diffusion 

In order to see that Yamazakura properly handles the resistive terms in the induction equation, we test a magnetic 
field diffusion problem in ID. Here, we solve not the full set of the MHD equations but the magnetic diffusion equation 
with a constant resistivity, 



dB _ d 2 B 

dt ^ dx 2 ' 



(ii) 



by freezing fluid in the induction equation. Setting the initial condition with an anti-parallel magnetic field, B = —Bo 
for x < and B = Bq for x > 0, the exact solution at time t is written as 



2Bn [ x l^ 1 



V 7 *" Jo 

B{x,t) = *% f 
Jo 



e u du for x < 0, 



x/ y/Arjt 



e u du for x > 0. 



(12) 



We set £?o = 100 and n = 10. The calculation is done with 800 numerical cells for x G [—5,5]. The Results are 
displayed in Fig. 61, in which one can find good agreement between the numerical and exact solutions. 
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Figure 61. Results of the magnetic field diffusion calculation. The numerical solutions are shown by the crosses while the analytic 
solutions are by the solid lines. 

Spontaneous Fast Reconnection 

We finally check that Yamazakura properly deals with the multi-dimensional resistive-MHD equations, by the sim- 
ulation of a spontaneous fast reconnection firstly done by Ugai (1999). The initial current system is constructed as 
follows: 

r sin(7n//2) for y < 1, 

1 for 1< y < 3.6, 

cos[(y - 3.6)7r/1.2] for 3.6 < y < 4.2, 
for 4.2 < y, 

l-B x (-y) /y/4nr for y < 0, 
p = (l+A,-^/47T)/2, 

p = 2p/(l + p ), 



B x (y)/y/*v = 



(13) 

where /?o = 0.15. A velocity and the x, y-component of magnetic field are initially zero. The adiabatic index is taken 
as 7 = 5/3. A reconnection of magnetic field- lines is initiated by an anomalous resistivity; 



r](r,t)- 



(k R [V d (r,t)-V c ] for V d >V c , 
\0 for V d < V c , 



where 



V d (r,t) = \J(r,t)/p(r,t)\, 
2p 



V c = V c n 



(14) 
(15) 

(16) 



The parameters are taken as = 0.03, V c ,o = 4, and a = 0.5. For details of the resistivity model, see Ugai (1999). 
The anomalous resistivity model (14) is assumed for t > 4. During < t < 4, a localized resistivity bellow is imposed 
to disturb the initial ID-configuration: 



r 1 (r) = r 1d exp[-{x/l.l) 2 - (y/1.1) 2 



(17) 



where rjd = 0.02. The computation is done with 800 x 800 numerical cells covering the domain of (x,y) G [—20,20] x 
[-6,6]. 

Some important results of the computation are shown in Fig. 62 and 63. The global magnetic-field distributions are 
displayed in Fig. 62 for t=18 and 30. We found that they are similar to those obtained by Ugai (1999), and also those 
of Feng et al. (2006). The left panel of Fig. 63 represents profiles of v x and B y along the x-axis at t = 24 and t = 30, 
while the right panel shows the evolution of the resistivity n and electric field E at the origin, v y at (x,y) = (0,0.9), 
and a magnetic flux $ evaluated by 



/ B x 

Jy>0 



(x = 0,y)dy. 



For these results, we also find a rough agreement with the above two works, although there are some insignificant 
differences. 

With all these results presented in the appendix, we expect that the present MHD collapse simulations done by 
Yamazakura offer proper results. 
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Figure 62. Magnetic field distributions at t=18 and 30 for the calculation of spontaneous fast reconnection. 




Akiyama, S., Wheeler, J. C, Meier, D. L., & Lichtenstadt, I. 2003, ApJ, 584, 954 

Auriere, M., Donati, J.-F., Konstantinova-Antova, R., et al. 2010, A&A, 516, L2 

Balbus, S. A. 1995, ApJ, 453, 380 

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 

Balsara, D. S., & Spicer, D. S. 1999, J. Chem. Phys., 149, 270 

Braithwaite, J. 2009, MNRAS, 397, 763 

Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 

Brio, M., & Wu, C. C. 1988, Journal of Computational Physics, 75, 400 

Burrows, A., Dessart, L., Livne, E., Ott, C. D., & Murphy, J. 2007, ApJ, 664, 416 

Diehl, R., Halloin, H., Kretschmer, K., et al. 2006, Nature, 439, 45 

Donati, J.-F., Howarth, I. D., Bouret, J.-C, et al. 2006, MNRAS, 365, L6 

Endeve, E., Cardall, C. Y., Budiardja, R. D., et al. 2012, ApJ, 751, 26 

Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 

Feng, X., Hu, Y., & Wei, F. 2006, Sol. Phys., 235, 235 

Ferrario, L. & Wickramasinghe, D. T. 2006, Mon. Not. Astron. Soc, 367, 1323 
Foglizzo, T., Scheck, L., & Janka, H.-T. 2006, ApJ, 652, 1436 

Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509 
Guilet, J., Foglizzo, T., & Fromang, S. 2011, ApJ, 729, 71 



Effects of Resistivity on Magnetized Core-Collapses 



Gustafsson, I. 1983, "Modified incomplete Cholesky (MIC) methods", in: D.J. Evans, ed., Preconditioning Methods; Theory and 

Applications (Gordon and Breach, New York), 265-293. 
Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 

Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ, 435, 339 
Iwakami, W., Kotake, K., Ohnishi, N., Yamada, S., & Sawada, K. 2009, ApJ, 700, 232 
Kotake, K., Sawai, H., Yamada, S., &; Sato, K. 2004, ApJ, 608, 391 
Kurganov, A., &; E. Tadmor 2000, Journal of Computational Physics, 160, 241 
Kuroda, T., & Umeda, H. 2010, ApJS, 191, 439 
LeBlanc, J. M., & Wilson, J. R. 1970, ApJ, 161, 541 
Leahy, D., & Ouyed, R. 2007, arXiv:0710.2114 
Liebendorfer, M. 2005, ApJ, 633, 1042 

Liska, R., & Wendroff, B. 2003, SIAM J. Sci. Comput. 25, 995 
Marek, A., Janka, H.-T., & Miiller, E. 2009, A&A, 496, 475 
Masada, Y., Sano, T., & Takabe, H. 2006, ApJ, 641, 447 

Moiseenko, S. G., Bisnovatyi-Kogan, G. S., & Ardeljan, N. V. 2006, MNRAS, 310, 501 

Murphy, J. W., & Burrows, A. 2008, ApJ, 688, 1159 

Nordhaus, J., Burrows, A., Almgren, A., & Bell, J. 2010, ApJ, 720, 694 

Obergaulinger, M., Aloy, M. A., & Miiller, E. 2006, A&A, 450, 1107 

Obergaulinger, M., Cerda-Duran, P., Miiller, E., & Aloy, M. A. 2009, A&A, 498, 241 

Obergaulinger, M., & Janka, H.-T. 2011, arXiv:1101.1198 

Papaliolios, C, Krasovska, M., Koechlin, L., Nisenson, P., &; Standley, C. 1989, Nature, 338, 565 
Sawai, H., Kotake, K. & Yamada, S. 2005, ApJ, 631, 446 
Sawai, H., Kotake, K. & Yamada, S. 2008, ApJ, 672, 465 

Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 1998a, Nuclear Physics A, 637, 435 

Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 1998b, Progress of Theoretical Physics, 100, 1013 

Shibata, M., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2006, Phys. Rev. D, 74, 104026 

Spitzer, L. 1956, Physics of Fully Ionized Gases, New York: Interscience Publishers 

Spruit, H. C. 2002, A&A, 381, 923 

Stone, J. M., & Norman M. L. 1992, ApJS, 80, 791 

Sumiyoshi, K., Yamada, S., & Suzuki, H. 2007, ApJ, 667, 382 

Sumiyoshi, K., Yamada, S., Suzuki, H., Shen, H., Chiba, S., & Toki, H. 2005, ApJ, 629, 922 
Symbalisty, E. M. D. 1984, ApJ, 285, 729 

Takiwaki, T., Kotake, K., Nagataki, S., & Sato, K. 2004, ApJ, 616, 1086 
Takiwaki, T., Kotake, K., & Sato, K. 2009, ApJ, 691, 1360 
Thompson, C, & Duncan, R. C. 1993, ApJ, 408, 194 
Thompson, T. A., Quataert, E., & Burrows, A. 2005, ApJ, 620, 861 
Ugai, M. 1999, Phys. Plasmas, 6, 1522 

Wade, G. A., Grunhut, J., Grafener, G., et al. 2012, MNRAS, 419, 2459 

Wang, L., Howell, D. A., Hoflich, P., & Wheeler, J. C. 2001, ApJ, 550, 1030 

Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914 

Yahil, A. 1983, ApJ, 265, 1047 

Yamada, S., & Sawai, H. 2004, ApJ, 608, 907 

Yoshizawa, A. 1990, Phys. Fluids B, 2, 1589 

Ziegler, U. 2004, Journal of Computational Physics, 196, 393 



